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I .  INTRODUCTION 


The  trend  in  packaging  of  electronic  components  is  to 
place  more  power  in  an  ever  decreasing  space.  Computers  that 
occupied  an  entire  room  20  years  ago,  now  occupy  a  small 
corner  of  that  room.  With  the  size  of  chips  decreasing,  the 
volumetric  power  generation  (q'")  is  of  the  order  of  MW/m^  in 
magnitude. 

This  can  lead  to  high  temperature  operation  of  the  chips 
with  a  resulting  degradation  in  performance  and  or  failure. 
The  military  has  set  the  maximum  temperature  of  a 
semiconductor  junction  to  between  100  and  110  'C  [Ref.  1]. 

This  investigation  deals  with  the  steady  state  temperature 
distribution  and  resulting  thermally  induced  stresses  in  a 
tri-material  electronic  package  as  shown  in  Figure  1.1, 
dimensions  not  shown  to  scale.  The  package  consists  of  a 
semi-conductor  (chip)  attached  to  a  substrate  material  by  a 
solder  joint.  The  chip  is  a  source  of  thermal  energy.  The 
package  is  exposed  to  a  convective  environment  on  all  sides. 


Figure  1 . 1  Tri-naterial  electronic  package 
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with  the  increase  in  volumetric  power  generation  the 
overall  temperature  of  the  package  is  increased.  The  stresses 
in  the  package  arise  from  two  sources,  (a)  thermal  strains 
(€^)  caused  by  temperature  gradients  within  a  material  and  (b) 
differences  in  the  materials  mechanical  properties, 
coefficient  of  thermal  expansion  (a)  and  stiffness  (E).  At 
material  interfaces,  displacement  continuity  requires 
expansion  or  contraction  accommodations  between  materials. 
There  are  two  such  interfaces  where  this  will  occur,  at  the 
chip/solder  interface,  and  at  the  solder/substrate  interface. 

A.  BACKGROUND 

Thermal  stresses  of  bi-material  assemblies  were 
investigated  by  Timoshenko  in  a  1925  paper  [Ref.  2],  where  a 
bi-metallic  thermostat  assembly  was  subjected  to  a  uniform 
temperature  field.  Timoshenko  obtained  analytical  expressions 
for  the  bending  stresses  through  the  cross  section  of  the 
bimetallic  strip.  He  also  noted  that  the  'distribution  of 
shearing  stresses  along  the  bearing  surface  cannot  be 
determined  in  an  elementary  way,  that  they  are  of  a  local  type 
concentrated  near  the  ends  of  the  strip ' ,  and  that  the 
shearing  stresses  can  be  of  the  same  magnitude  of  the  bending 
stresses.  He  mentions,  in  passing,  the  existence  of  'local' 
normal  stresses  between  the  assembly  interfaces  but  says 
nothing  of  their  magnitude.  Goland  and  Reissner  [Ref.  3] 
determined  the  shear  and  peeling  stresses  at  the  interface 
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of  a  cemented  lap  joint,  where  the  peeling  stress  tends  to 
pull  the  materials  apart. 

In  a  comprehensive  review  of  the  subject  of  tri-material 
assemblies,  Suhir  [Ref.  4]  took  the  analysis  one  step  further 
by  investigating  the  stresses  in  a  tri-material  system  due  to 
a  thermal  environment.  His  model  is  based  on  an  assembly 
fabricated  at  elevated  temperatures  and  subsequently  cooled. 
Specific  restrictions  were  placed  on  the  tri-material  assembly 
in  a  uniform  temperature  field.  He  obtained  analytical 
solutions  of  the  problem.  For  the  present  study  the  general 
tri-material  problem  is  solved  numerically.  In  particular  the 
tri-material  system  is  that  of  an  electronic  package. 

B.  PROBLEM  DEPIMITIOH 

The  problem  was  partitioned  into  two  parts,  a  thermal 
analysis,  whereby  the  heat  conduction  equation  was  solved,  and 
a  thermal  stress  analysis  whereby  the  stresses  produced  by  the 
thermal  field  were  obtained.  Both  problems  were  solved 
numerically  by  the  Finite  Element  Method  (FEM).  For  the 
numerical  analysis,  symmetry  is  invoked  along  the  centerline 
oi  the  system  as  shown  in  Figure  1.1.  The  steady  state 
thermal  study  focuses  on  the  effects  material  properties, 
geometric  configuration,  convective  cooling,  and  encapsulation 
have  on  temperature  distribution  in  an  electronic  package. 
Once  the  temperature  field  is  known  the  data  is  input  into  a 
program  that  assigns  temperatures  to  another  FEM  mesh  which  is 
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used  in  the  stress  analysis.  The  stresses  examined  include 
bending  stress  at  the  centerline  of  the  package  the 

bending  (<7,,^),  normal  (cr^)/  and  shear  (r^),  stresses  along  the 
chip/solder  interface,  and  the  bending  normal 

shear  (r^)  stresses  at  the  solder/substrate  interface. 
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II.  PEM  DEVELOPMENT  OP  THERMAL  PROBLEM 


A.  PROBLEM  STATEMENT  AND  ASSUMPTIONS 

As  noted  in  the  introduction/  the  purposes  of  this 
investigation  are  to  (a)  determine  the  variation  of 
temperature  in  the  electronic  package  and  (b)  determine  the 
stresses  associated  with  these  temperatures.  The  results  here 
are  the  steady  state  results . 

The  steady  state  thermal  study  focuses  on  effects  that 
certain  key  parameters  have  on  system  behavior.  They  include, 
the  effects  total  power  (P)  in  the  chip,  material 

combinations,  convective  cooling,  geometric  configuration  and 
encapsulation  have  on  the  temperature  distribution  in  the 
electronic  system. 

Finite  element  codes  were  used  to  solve  for  the 
temperature  distribution  and  stresses  of  the  electronic 
package.  In  order  to  limit  the  mathematical  complexity  of  the 
problem  the  following  approximations  and  specific  assumptions 
for  FEM  program  HEATSTEADY  in  Appendix  A.  are  made  as  follows: 

•  The  thermal  contact  resistances  due  to  'imperfect'  contact 
at  the  solder  interfaces  are  negligible. 

•  The  thermal  conductivity  and  coefficient  of  thermal 
expansion  (k  and  a)  are  considered  constant  over  the  range 
of  temperatures  encountered  in  the  study. 
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•  Because  the  temperatures  are  not  very  large  radiation  heat 
transfer  between  the  package  and  the  surrounding  enclosure 
is  negligible. 

•  The  convection  heat  transfer  coefficient  (h)  is  constant 
over  the  entire  edge  upon  which  it  acts. 

All  of  these  assumptions  are  reasonable  within  the  scope 
of  this  study.  Once  the  temperature  distribution  of  the 
package  is  obtained  the  position  and  temperature  of  all  the 
Global  nodal  points  are  entered  into  program  Filter  in 
Appendix  B.  to  obtain  temperatures  that  correspond  to  FEM 
program  Weld  in  Appendix  C.  to  determine  the  stresses 
developed  in  the  package.  The  development  of  program  WELD  is 
discussed  in  Chapter  III. 

B.  FEM  FORMULATION  OF  THERMAL  PROBLEM 

The  Galerkin  FEM  is  an  approximation  method  which 
transforms  a  linear  partial  differential  equation  into  a 
system  of  linear  algebraic  equations.  Using  the  two- 
dimensional  heat  equation  [Ref.  5]: 

V  •  (kV  T)  +  q'"  =  0  (2.1) 


and  boundary  conditions 

Cauchy:  -k-|^  =  h(T-TJ 

an 


(2.2) 


Nuemann : 


(2.3) 
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where  the  V  operator  denotes 


i—  +  7-^ 

dx  ■'  dy 


,  and  thus 


V-(JtVT)  =  icV^T  =  4—  +  — 

[dx^  dy^ 

where  k  is  the  thermal  conductivity,  T  is  the  temperature,  and 
q'"  is  the  volumetric  heat  generation.  The  temperature 
distribution  can  be  determined  at  system  nodal  points  of  the 
electronic  package.  For  this  method,  a  linear  triangular 
shape  function  (N,-)  which  possesses  the  Kronecker  Delta 
property  (equation  2.5); 

N^inodej)  =  if  2  =  j  then  =  1 

if  i  *  j  then  -  0  '  * 

is  used.  The  linear  shape  functions  maintain  function 

continuity  throughout  the  domain. 

Heat  flux  continuity  at  material  interfaces  is  built  into 
the  FEM  formulation.  The  thermal  FEM  grid  and  boundary 
conditions  are  given  in  Figure  2.1  and  Figure  2.2 
respectively . 
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B  DRY  BOUNDARY  CONDITION 


D 

-k  4^  =  h  (T-T_  ) 
d  y  00 

B 

-k  4^  =  h  (T-T^) 
ox  00 

C 

-k  4^  =  h  ) 

d  y  CD 

B 

E 

B 

ox 

BOUNDARY  CONDITIONS  FOR 
THERMAL  PROBLEM 


Figure  2 . 2  FEM  thermal  problem  boundary  conditions 


An  approximate  solution,  t,  for  temperature,  T(x,y),  is 
formed  as  follows: 

T  »  B  =  {N\'^[f}  (2.6) 

where  T  is  the  exact  solution  of  the  heat  equation  in 
continuous  space,  t  is  the  approximate  solution  in  discrete 
space,  {N}^  is  the  transpose  of  a  column  vector  of  the 
triangular  linear  shape  functions,  and  {T}  is  the  vector  of 
nodal  temperatures . 

After  the  approximation  is  formulated,  the  next  step  is  to 
form  the  Residual,  R,  as  follows: 

R  =  ^iC)  -  g"'(s)  (2.7) 

where  q'"  is  the  heat  generation  term  and  SP  denotes  the 
differential  operator  which  in  the  case  of  the  heat  equation 
is  defined  by: 

S£(t)  =V‘{kVt)  (2.8) 

With  this  substitution,  the  residual  becomes: 

R=  -  q'"[s)  (2.9) 

From  the  residual,  the  Galerkin  Equations  are  formed: 

r  {N){R)d3  =  {0)  (2.10) 

J  D 
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where  D  denotes  surface  integration  and  where  {0}  is  the  null 
vector.  Further  substitution  for  R  into  the  Galerkin  vector 
equation  results  in: 

f  {N)[V'{k^{{N)^{T)))  ds  -  f  {N}q'"(3)ds  =  {o}  (2.11) 

J  V  J  D 

To  solve  the  Galerkin  Equation/  Green's  Theorem  is  evoked 
which  yields: 

(f  {n)  (kVT)  dB  -  f  [V{n}  (kV({N}^{T}) )  1  ds  +  f  {N}(j'"ds=  {o} 

is  J  D  J  D 

where  B  subscript  on  the  integral  denotes  evaluation  of  these 
integrals  around  the  boundary  of  domain  D  of  the  electronic 
package.  When  integrated  the  boundary  integral  and  the  heat 
generation  integral  form  the  “force"  vector  {F}.  The  middle 
integral  forms  the  [A]  matrix.  Thus  the  heat  conduction 
differential  equation  becomes/ 

[A]  (r)  =  (F)  (2.13) 

where  {T}  contains  the  vector  of  nodal  temperatures. 

C.  THERMAL  PROGRAM  DESCRIPTION  AMD  VALIDATION 

The  main  thermal  program/  HEAT. FOR/  a  VAX  Fortran  77  Code 
for  the  FEM  thermal  analysis  is  contained  in  Appendix  A.  The 
program  begins  by  reading  in  a  data  file  DATA3.DAT.  This 
input  data  is  then  processed  through  several  subroutines  that 
perform  the  FEM  analysis. 
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1.  Subroutine  Input  Description 

The  following  information  is  supplied  to  the  Input 

subroutine  from  a  data  file. 

Input  file  form:  NC0L1,NC0L2/NR0W, 

XPOS , ypos , DIVX, DIVY, 

THERCON , GEN , QTPR , 

DIVSUB, DIVXAIR, DIVYAIR, 

NBCLOW, NBCUP , NBCVERT, NBCSUB, TAMB , 

BLOW , HUP , HVERT , HSUB 


Parameter 

NCOLl 

NCOL2 

NROH 

XPOS 

YPOS 

DIVX 

DIVT 

THERCON 

GEN 

QTPR 

DIVSUB 

DIVXAIR 

DIVYAIR 

NBCLOW 

NBCUP 


Description 

The  number  of  columns  of  substrate  in  the  FEM  grid. 

The  number  of  columns  of  chip  in  the  FEM  grid. 

The  number  of  rows  in  the  FEM  grid. 

Starting  position  for  new  grid  density  along  X 
axis . 

Starting  position  for  new  grid  density  along  Y 
axis . 

Spacing  of  grid  density  along  X  axis. 

Spacing  of  grid  density  along  Y  axis. 

Thermal  conductivity  (k). 

0  =  no  heat  generation  term 
1  =  heat  generation  term 

Amount  of  volumetric  heat  generation  ((^")» 

Number  of  rows  in  FEM  grid  composed  of  substrate. 

Number  of  columns  in  FEM  grid  composed  of  air. 

Number  of  rows  in  FEM  grid  composed  of  air. 

0  =  insulated  boundary  at  the  bottom  of  substrate. 
1  =  convective  boundary  at  the  bottom  of  substrate. 

0  =  insulated  boundary  at  the  top  of  chip. 

1  =  convective  boundary  at  the  top  of  chip. 
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NBCVERT 


0  =  insulated  boundary  on  the  side  of  the  chip. 

1  =  convective  boundary  on  the  side  of  the  chip. 


NBCSUB  0  =  insulated  boundary  on  top  of  the  substrate. 

1  =  convective  boundary  on  top  of  the  substrate. 

HLOW  Convection  coefficient  corresponding  to  NBCLOW. 

HUP  Convection  coefficient  corresponding  to  NBCUP. 

HVERT  Convection  coefficient  corresponding  to  NBCVERT. 

HSUB  Convection  coefficient  corresponding  to  NBCSUB. 

2 .  Subroutine  Grid  Description 

Subroutine  GRID  takes  the  input  information  and 
constructs  a  grid  of  horizontal  and  vertical  lines  on  a 
geometry  as  shown  in  Figure  2.1.  This  grid  will  be  used  to 
form  the  triangular  elements  used  in  the  thermal  FEM  analysis. 
This  subroutine  allows  the  user  to  generate  a  variety  of 
meshes  and/or  refinements  and  permits  validation  of  grid 
independence  very  quickly. 

To  optimize  computer  time  the  subroutine  is  able  to 
generate  different  grid  meshes  with  a  variety  of  spacing 
options  as  shown  in  Figure  2.3.  Fine  meshes  are  used  at 
convective  boundaries  and  material  interfaces.  A  coarse  mesh 
is  used  within  a  material  away  from  boundaries.  The 
subroutine  allows  the  user  to  independently  vary  the  mesh 
density  both  horizontally  and  vertically  for  up  to  nine 
different  densities  in  either  direction.  The  number  of  mesh 
densities  can  be  increased  by  simply  adding  a  number  of  "IF 
THEN"  statements  at  the  top  of  the  program  (line  34). 
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fine  coarse  coarse 


fine 

Figure  2.3  Fine  and  coarse  vertical  mesh  and  fine  and  coarse 
horizontal  mesh 

Local  nodal  points  are  assigned  to  each  element,  one 
at  each  corner,  counterclockwise  as  shown  in  Figure  2.4. 


Figure  2 . 4  Triangular  element  numbering 

Next  the  subroutine  generates  a  correspondence  table 
which  identifies  the  correspondence  between  local  nodal  points 
and  global  nodal  points.  Next,  each  triangular  element  is 
identified  with  a  specific  material  and  assigned  a  thermal 
conductivity  (k)  and  whether  or  not  it  will  have  a  volumetric 
heat  generation  term  (q'").  Lastly  boundaries  are  assigned 
according  to  Figure  2.2. 
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3.  Subroutin*  H««taat,  BMtmod  and  Lsarg  Dascription 

Subroutine  HEATMAT  constructs  the  global  matrix  [A] 
and  the  heat  generation  vector  {F}  according  to  the 
correspondence  table  made  in  subroutine  GRID. 

Subroutine  HEATMOD  is  not  used  in  this  program  but  can  be 
used  to  modify  the  global  matrix  [A]  for  a  Dirichlet  boundary. 
HEATMOD  replaces  the  Galerkin  equation  for  the  specific  nodal 
points  with  the  specified  temperatures.  LSARG  is  an  equation 
solver  provided  by  the  IMSL  Library  used  to  solve  equation 
2.13. 

4.  Subroutine  Output  Description 

Subroutine  OUTPUT  creates  data  files  and  performs  an 
energy  balance  between  the  thermal  energy  which  is  convected 
at  the  boundaries  and  the  thermal  energy  generated  by  the  chip 
according  to: 

^  h,  {T,  -  TJL,  = 


where  h^  and  1^  correspond  to  the  convection  coefficient  and 
length  of  each  convective  boundary,  T^  corresponds  to  the 
average  temperature  across  the  boundary  and  Ag^^p  is  the  area 
of  the  chip.  This  heat  balance  is  conducted  for  each  case 
study. 
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D.  VALIDATION  I.  DIRICHLET  AND  CAUCHY  BOUNDARY 

A  two  dimensional  body  subjected  to  both  Dirichlet  and 
Cauchy  boundary  conditions  was  taken  from  [Ref.  6]  and  run  on 
program  HEAT  to  validate  the  programs  ability  to  handle  both 
boundary  conditions.  The  results  of  program  HEAT  were 
identical  to  Reference  6. 

E.  VALIDATION  II.  QRID  INDEPENDENCE 

In  order  to  ensure  that  the  numerical  errors  are  kept  to 
a  minimum  the  size  of  the  mesh  was  decreased  with  incremental 
steps  until  there  was  a  less  than  1%  change  in  the  solution 
from  one  grid  to  the  next  finest  grid  (grid  independence) . 

The  final  mesh  chosen  has  125  degrees  of  freedom  (DOF) 
with  202  elements.  The  mesh  which  provided  less  than  1% 
difference  had  161  DOF  with  268  elements.  The  total  power  of 
the  chip  used  is  P  ==  40  W.  The  total  power  convected  for  the 
DOF  =  125  grid  was  P<,onv  =  39.95  W.  The  total  power  convected 
for  the  DOF  =  161  grid  was  P^^^^  =  39.94  W.  The  percent 
difference  in  total  power  P  from  the  coarse  mesh  to  the  fine 
mesh  is  0.03%. 

F.  VALIDATION  III.  ENERGY  BALANCE 

As  mentioned  previously,  there  is  an  energy  balance  in 
subroutine  OUTPUT  which  is  given  by  equation  2.14.  This 
checks  program  HEAT'S  ability  to  correctly  handle  a  generation 
term.  The  total  power  for  the  majority  of  case  studies  was 
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equal  to  P  =  1  W.  The  average  convective  power  output  was 

q  =  .99785  W  for  a  percent  difference  of  0.21%. 

^conv 
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III.  FEM  STRESS  DEVELOPMENT 


The  finite  element  stress  formulation  utilizes  a  recently 
developed  element  which  provides  for  axial  and  lateral 
displacement  continuity.  The  results  of  the  stress  code  are  in 
good  agreement  with  existing  solutions  [Ref.  2,4]  to  problems 
with  uniform  temperature  fields. 

A.  FORMULATION  OF  THE  MODEL 

A  brief  description  of  the  FEM  formulation  for  stresses 
follow.  In  Figure  3.1,  each  element  has  six  degrees  of 
freedom;  axial  displacements  at  the  four  corner  points,  and 
lateral  displacements  at  the  two  ends.  An  advantage  of  the 
element  is  that  axial  and  lateral  displacement  continuity 
results . 


Figure  3.1  A  typical  element  with  6  degrees  of  freedom 
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The  axial  displacement  field  a(x,y)  is  assumed  to  be  linear  in 
both  the  axial  and  transverse  directions.  That  is, 

u(x,y)  =  '^Ni(x)[H^{y)  Ui  +H2(y)  Ui]  (3.1) 


Superscripts  b  and  t  on  the  nodal  displacements  refer  to 
bottom  and  top  displacements  respectively.  The  linear  shape 
functions  N-  and  are  given  by 

Hi  W  ‘  ^  Hj  =  ^  N,  (X)  =  W,  = 

The  lateral  displacement  field  is  given  by 

vix)  =^Ni(x)v^  (3.2) 


The  strain-displacement  relations  are, 


'  ^  -  Ul'l 


(3.3) 


and 


rxy 


=  ^[“l'  -  “ll  *  -^[“2'  -  ^2“]  * 


Vj  -  V, 


(3.4) 
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Defining  the  axial  displacement  vector  as 


=  Ui"  ui)  (3.5) 

and  the  stiffness  matrix  as 

[K^]  =  (3.6) 


where  the  B  vector  is, 


= 


dx 


dx 


dx 


(3.7) 


and  the  force  vector  due  to  temperature  as. 


(F  }  s  ^  f^f^lBlEaATdydx 

olJQ  Jo 


(3.8) 


gives  the  bending  matrix  equations  as, 

(3.9) 

Equation  (3.9)  defines  the  bending  behavior.  Behavior  due  to 
shear  is  obtained  as  follows.  Define  the  row  vector  of 
displacement  degrees  of  freedom  as 

=  V,  ui  u/  v^)  (3.10) 

The  shear  stiffness  matrix  is  given  by 
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[Kg]  =  J^^f^^{B')G{B')dydx 


(3.11) 


where 


(B'>  = 

dy 


ay 


ajy^ 

ax 


a/fi 

ay 


AT, 


a//,  ajv; 


_  2)  (3.12) 

ay  ax 


which  gives  the  equations  for  shear  behavior 

[JCj{fi^}  =  {0}  (3.13) 

The  matrix  equations  for  bending  and  shear  behavior  are 
combined  to  give  the  stiffness  equations  for  the  system. 

B.  PRELIMIBARl  RESULTS 

Assessment  of  the  present  model  was  obtained  from 
comparisons  with  previous  analytical  models  [Ref.  2,4].  Exact 
agreement  for  bending  stresses  was  obtained  when  the  present 
model  was  applied  to  a  sample  bimetallic  case  solved  in 
[Ref.2].  For  a  tri-metallic  case  the  present  model  produced 
shear  stress  results  in  reasonably  good  agreement  with  results 
obtained  from  [Ref.  4]. 
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IV.  STEADY  STATE  TEMPERATURE  AMALYSIS 


A.  OVERVIEW  OF  STEADY  STATE  TEMPERATURE  ANALYSIS 

The  analysis  begins  with  an  FEM  approximation  of  the 
temperature  distribution  of  the  tri-material  electronic 
package.  First  the  chip  volumetric  power  generation  is 
examined  to  determine  the  effect  an  increase  in  q'"  has  on  the 
chip  temperature.  Keeping  the  power  of  the  chip  and  the  level 
of  convective  cooling  constant,  several  packaging  materials 
are  combined  to  determine  how  the  temperature  distribution 
varies  with  the  variation  in  the  material  thermal 
conductivities.  This  is  followed  by  a  convection  cooling 
analysis,  in  which  temperatures  and  temperature  gradients  are 
examined  for  a  common  electronic  package  while  the  convection 
coefficient  h  varies.  The  package  geometry  is  then  studied, 
holding  the  overall  package  height  equal  and  varying  the  chip 
size  as  a  percentage  of  overall  package  height.  This  is  done 
for  the  case  of  constant  heat  flux,  and  for  the  case  of 
constant  volumetric  heat  generation.  Lastly  the  effects  of 
placing  a  protective  coating  "encapsulation"  on  the  chip  is 
studied  to  determine  the  variation  in  temperature  distribution 
from  a  chip  which  has  not  been  encapsulated. 
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B.  EFFECT  OF  VOLUMETRIC  POWER  GEMERATIOM  OR  TEMPERATURE  OF 


THE  CHIP 

The  heat  fluxes  of  electronic  packages  is  expected  to 
reach  nearly  10^  W/m^  by  the  year  2000.  This  is  two  magnitudes 
in  power  over  the  chip  studied  here.  The  total  power  output 
(P)  of  the  chip  studied  is  one  watt.  This  translates  to  a 
volumetric  power  generation  (q'")  of  20  MW/m^  for  a  device  5mm 
X  5mm  X  2mm  in  dimension  or  a  heat  flux  (q")  of  4x10*  W/m^. 
The  dimensions  chosen  are  of  a  typical  electronic  device  size. 

The  temperature  results  throughout  the  study  can  be 
modified  for  any  power  output  due  to  the  linearity  of  the 
generation  term  in  the  heat  equation.  In  Figure  4.1  the 
temperature  of  the  chip  is  plotted  against  chip  volumetric 
power  output  (q''').  All  other  quantities  must  be  constant 
for  this  to  remain  true,  including  boundary  conditions  and 
other  materials  within  the  package.  For  the  particular 
geometry  studied,  temperature  (T)  as  a  function  of  volumetric 
power  output  q'"  W/m^  is  given  by  equation  4.1 

T  =  m  *  q"'  +  r.  (4.1) 

where  the  change  in  temperature  per  unit  change  in  volumetric 
power  output  is  b  =  4.65  *C/MW/m^  and  T,  is  the  ambient 
temperature  of  the  convective  fluid. 
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0  5  10  15  20 

VOLUMETRIC  POWER  (q'"  MW/m*) 


Figure  4 . 1  Chip  temperature  as  a  function  of  rolumetric  power 
C.  EFFECT  OF  VARIOUS  MATERIAL  COMBINATIONS 

The  materials  used  in  this  investigation  along  with  their 
key  thermal  properties  are  listed  in  Table  4.1.  The  set  of 
materials  that  comprise  the  electronic  tri-material  package 
focused  on  for  the  majority  of  this  study  are  the  silicon 
chip,  a  Pb-Sn  solder,  along  with  several  substrate  materials 
including.  Epoxy  Fiberglass,  Polyimide  Fiberglass,  Alumina, 
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and  Aluminum  Nitride.  These  are  the  most  common  materials 
used  in  electronic  packaging.  The  other  chip  and  solder 


materials 

considered  in  this 

study  account  for  large 

variations 

in  thermal  conductivity  and/or 

thermal  expansion 

coefficient  (a). 

TABLE  4.1 

ELECTRONIC  PACKAGE  MATERIALS  AMO 

PROPERTIES 

Use 

material 

a  (10-^/‘C\ 

Chip 

Silicon 

150 

27 

Gallium  Arsenide 

58 

57 

Substrate 

Epoxy  Fiberglass 

0.16 

140  -  180 

Polyimide  Fiberglass 

0.35 

120  -  160 

Aliimina 

18 

60 

Aluminum  Nitride 

230 

33 

Solder 

Pb(5%)/Sn 

63 

290 

Au(20%)/Sn 

57 

159 

Air  Gap 

0.026 

- 

Throughout  the  study  within  this  section  the  following 
parameters  are  held  constant;  volumetric  power  (q'")  of  the 
chip,  ambient  temperature  T,  of  the  convective  fluid,  and 
coefficient  of  convection  cooling  (h) .  Figure  4.2  shows  the 
non-dimensional  temperature  profiles  of  the  package  at  a  non- 
dimensional  distance,  =  0.564.  The  plotted  values  of  T^ 
were  obtained  at  the  FEM  grid  nodes.  The  temperatures  which 
lie  along  =  .564  were  chosen  because  they  pass  through  the 
solder  joint  of  the  package.  Thermal  stresses  will  result  not 
only  from  temperature  gradients  within  a  material,  but  also 
from  the  direct  contact  of  the  three  different  materials 
within  the  package  which  have  large  differences  in  their 
coefficients  of  thermal  expansion  (a).  In  Figure  4.2  the 
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abscissa  is  a  normalized  temperature  given  by  equation  4.2, 


T  = 


(4.2) 


where  is  the  FEM  grid  temperature.  The  ordinate  is  the 
non-dimensional  y  location  (py  that  is  given  by  equation  4.3. 

where  H  is  the  total  height  of  the  package  and  y  is  the  FEM 
grid  global  nodal  location. 


NON  DIMENSIONAL  TEMPERATURE  T^= 

CO 


Figure  4.2  Mon-diaenslonal  teaperature  for  alllcon  chip,  Pb-Sa 
solder  and  various  substrates. 
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All  sixteen  combinations  of  materials  were  analyzed  for  a 
temperature  profile.  There  was  no  significant  difference  in 
temperature  profile  with  regard  to  which  chip  or  solder 
material  was  used  in  the  package.  For  the  volumetric  power 
level  used  in  this  analysis  the  difference  in  thermal 
conductivity  of  the  Silicon  chip  (k*150  W/m'C)  or  Ga/As  chip 
(k=58  W/m‘C)  does  not  influence  the  temperature  profile.  A 
study  was  conducted  in  which  the  volumetric  power  was  raised 
by  four  times  the  initial  amount.  This  failed  to  produce  any 
significant  difference  in  temperature  profile  or  magnitude. 

The  temperature  profiles  show  that  the  substrate  is  the 
governing  material  which  determines  the  temperature  profile 
a  given  electronic  package  takes  on.  The  following  symbology 
will  be  used  in  further  analysis.  A  subscript  of  1  refers  to 
the  chip,  subscript  2  refers  to  the  solder,  and  subscript  3 
refers  to  the  substrate.  Low  thermal  conductivity  substrates 
used  in  this  analysis  have  ratios  of  kj/k^  <  0.0025.  High 
thermal  conductivity  substrates  used  in  this  analysis  have  a 
ratio  range  of  0.10  <  kj/k^  <  1.6. 

A  characteristic  of  a  high  kj/k^  package  is  a  uniform 
temperature  field  throughout  all  materials  in  the  package. 
All  materials  within  the  package  can  conduct  heat  very  well. 
With  the  low  kj/k^  substrates  a  temperature  gradient  is 
developed  within  the  substrate  in  both  the  X  and  Y  directions. 
As  a  result  of  not  being  able  to  conduct  the  heat  effectively 
through  the  low  conductive  substrate,  the  chip  temperature  is 
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Increased  by  approximately  40%  over  a  package  with  a  highly 
conductive  substrate. 


The  solders  have  thermal  conductivities  which  are  close  to 
the  chip  thermal  conductivity  (kj/k,  =  0.446)  and  do  not  effect 
the  temperature  profile.  However,  the  coefficient  of  thermal 
expansion  a  between  the  two  solders  chosen  vary  by  55%  thereby 
influencing  the  thermal  strain  profile  significantly.  The 
solder  joint  takes  on  the  chip  temperature  in  all  cases 
because  of  its  ability  to  conduct  the  heat  generated  by  the 
chip. 

Thermal  strain  profiles  are  presented  in  Figures  4.3  and 
4.4  for  the  two  different  solders.  In  these  figures,  the 


THERMAL  STRAIN  €_»  (  T-T  )  aCxIo’ni/m) 

I  CD 

Figure  4.3  Themal  strain  profile  for  silicon  chip,  Pb-8n 
solder  and  various  substrates. 
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NON  DIMENSIONAL  Y  LOCATION 


abscissa  is  the  difference  of  the  nodal  temperature  and 
initial  temperature  multiplied  by  the  coefficient  of  thermal 
expansion  a  given  by 

€r  =  {T,  -  TJ  *  a  (4.4) 

and  the  ordinate  is  given  by  equation  4.3. 

1.0 

m 

.8 

.6 


.4 

.2 

0 

0  .4  .8  1.2  1.6  2.0 

THERMAL  STRAIN  *  (T-T^)  a  (xlO*  m/m) 

T  ® 

Figure  4.4  Thermal  strain  profile  for  silicon  chip,  Au-Sn 
solder  and  various  substrates. 


1.  Low  Thermal  Conductivity  Substrate  Packages 


Further  analysis  is  conducted  by  dividing  the  packages 
into  two  groups,  low  thermal  conductivity  substrates  and  high 
thermal  conductivity  substrates.  The  two  substrates  chosen. 
Epoxy  Fiberglass  (kj/k^  =  0.00106)  and  Polyimide  Fiberglass 
(kj/k^  =  0.0023),  are  part  of  a  group  of  substrates  that  have 
kj  ranging  from  0.12  W/m*C  to  0.35  W/m*C  [Ref.  1].  The 
temperature  achieved  throughout  the  chip  for  each  of  the  case 
studies.  Figure  4.2,  are  approximately  constant,  that  is,  not 
dependent  on  position.  The  variation  in  temperature  of  the 
chip  was  limited  to  less  than  1.2%  over  the  entire  cross 
section.  The  solder  achieved  the  same  temperature  as  the  chip 
due  to  the  solders  high  thermal  conductivity.  The  variation 
between  the  chip  temperature  and  the  solder  temperature  is 
less  then  0.1%. 

Within  the  substrate  temperature  gradients  0^  and  0^ 
are  developed  in  both  the  X  and  Y  directions  respectively. 
Temperature  gradients  are  given  by; 


dT  _  0  _  (Tj  ~ 
dx  *  ”  (Xi  -  x^.i) 


(4.5) 


BT  _  0  _  (Yj  ~ 

dy  ~  ^  ~  iVj  -  yj*i) 


(4.6) 


where  T^  and  T^  are  FEM  Global  nodal  temperatures,  x^  and  yj 
are  FEM  Global  nodal  distances  in  the  X  and  Y  directions 
respectively. 
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These  temperature  gradients  are  only  found  in  low 
thermal  conductivity  substrates  and  an  explanation  for  this 
can  be  found  by  examining  Fourier's  Law,  equation  4.7, 

(4.7) 

where  q"  is  the  heat  flux  in  the  X  direction. 

The  solder  acts  as  a  sort  of  'heat  funnel'  for  the 
chip  to  transfer  its  thermal  energy  to  the  substrate  since  the 
air  pocket  under  the  chip  acts  as  an  insulator.  By  examining 
Fourier's  law  there  are  two  ways  to  increase  the  heat  flux  q". 

First  by  increasing  the  thermal  conductivity  (k)  of 
the  material,  and  secondly  by  increasing  the  temperature 
difference  (or  temperature  potential).  The  convection  heat 
transfer  from  the  chip  surface  to  the  surrounding  medium  is 
unable  to  transfer  all  the  thermal  energy  in  the  chip 
generated  by  the  power  source  q'"  and  therefore  results  in  the 
conduction  of  heat  flux  through  the  solder.  When  the  thermal 
conductivity  of  the  substrate  is  low  (as  in  the  Fiberglass 
substrates ) ,  the  temperature  difference  through  the  substrate 
must  become  proportionally  larger  to  account  for  the  heat  flux 
out  of  the  chip.  For  this  reason,  when  the  substrate  has  a 
high  thermal  conductivity,  the  temperature  potential  need  only 
be  very  small  to  conduct  the  same  amount  of  thermal  energy. 

In  the  X  direction  the  gradient  is  largest  at  the 
upper  right  hand  side  of  the  solder/substrate  interface.  As 
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Figure  4.5  shows,  the  gradient  6^  diminishes  by  roughly  72% 
when  only  25%  into  the  substrate.  The  temperature  gradient 


Figure  4.5  9,  in  epoxy  fiberglass  substrate 

has  decreased  by  88%  when  half  way  through  the  substrate. 
This  behavior  can  be  accounted  for  by  the  fact  the  thermal 
resistance  of  the  substrate  in  conduction  is  greater  then  the 
thermal  resistance  of  the  convection  surface  at  the  top  of  the 
substrate.  Equivalent  thermal  resistances  in  conduction  and 
convection  are  given  by; 
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R 


cond 


=:  Ai[ 

k  A 


(4.8) 


where  Ax  is  the  distance  between  two  points  of  interest 


conv 


1 

h  A 


(4.9) 


and  A  is  the  area  in  which  the  heat  flux  is  acting. 

Taking  the  area  under  the  solde**  joint  for  conduction 
and  the  top  of  the  substrate  to  the  right  of  the  solder  joint 
for  convection  the  following  equivalent  thermal  resistances 
are  obtained: 


,  ^  _ 0.002 _ 

0.16  (0.00088)  (1) 


14.2  “C/fiT 


(4.10) 


R 


conv 


_ 1 _ 

200  (0.001)  (1) 


5  ^C/W 


(4.11) 


Which  shows  clearly  that  the  thermal  energy  has  approximately 
1/3  less  resistance  at  the  convective  boundary  versus 
conduction  through  the  substrate. 

The  gradient  in  the  X  direction  on  the  left  hand  side 
of  the  solder  joint  is  only  41%  of  that  on  the  right,  as  shown 
in  Figure  4.6.  This  again  is  a  result  of  a  higher  thermal 
resistance  in  conduction  within  the  substrate.  With  the 
insulated  boundary  condition  imposed  in  the  FEM  formulation, 
the  gradient  in  the  X  direction  at  X=0  and  X=L  is  necessarily 


zero. 


Y  LOCATION  THROUGH  SUBSTRATE  (mm) 


TOP  OF  SUBSTRATE 


Figure  4.6  in  epoxy  fiberglass  substrate 

The  temperature  gradient  in  the  Y  direction,  0^,  is 
largest  at  49.02  ’C/mm  at  the  solder/substrate  interface  as 
shown  in  Figure  4.7.  Away  from  this  interface  the  gradient 
decreases  to  approximately  20  *C/mm  and  remains  constant 
throughout  the  entire  substrate. 


Y  LOCATION  THROUGH  SUBSTRATE  (mm) 


Figure  4.7  9y  in  epoxy  fiberglass  substrate 


The  position  where  the  largest  temperature  gradients 
in  both  the  X  and  Y  directions  are  found  are  shown  in  Figure 
4.8.  These  are  the  points  where  maximum  heat  fluxes  are 
conducted  through  the  solder/substrate  interface. 
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Figure  4,8  Naxinun  tenperature  gradients  in  epoxy  fiberglass 
substrate 

2 .  High  Thermal  Conductivity  Substrate  Packages 

The  two  substrates  chosen  in  this  range  are  Alumina  k, 
=  18  W/m*C  (kj/k^  =  0.12)  and  Aluminum  Nitride  (AlN)  kj  =  230 
W/m*C  (kj/k,  =  1.53).  Referring  back  to  Figure  4.2  the  non- 
dimensional  temperature  distribution  along  =  .564  for 

Alumina  and  AIN  graph  is  approximately  a  vertical  straight 
line.  Therefore  there  are  no  temperature  gradients  anywhere 
within  the  package.  The  reason  for  no  temperature  gradients 
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in  the  substrate  is  due  to  the  substrate's  high  thermal 
conductivity,  enabling  it  to  diffuse  a  large  heat  flux  without 
developing  a  temperature  gradient,  as  shown  by  Fourier's  Law 
of  heat  conduction,  eguation  4.7.  Figure  4.3  shows  the 
largest  shear  strains  at  the  solder  interface  have  been 
reduced  by  40%  over  the  low  thermal  conductivity  substrates 
due  to  the  increased  ability  of  the  high  thermal  conductivity 
substrates  to  conduct  heat  away  from  the  chip  thereby  lowering 
the  overall  temperature  of  the  package  and  decreasing  the 
thermal  expansion. 

There  is  however  no  appreciable  temperature  difference 
between  the  Alumina  and  AlN  substrate  despite  the  thermal 
conductivity  of  the  Aluminum  Nitride  being  approximately  13 
times  as  large  as  the  Alumina  thermal  conductivity.  The 
limiting  factor  for  reduction  of  temperature  of  the  package 
has  become  the  convective  thermal  resistance,  which  must  be 
lowered  in  order  to  lower  the  package  temperature.  The 
conductive  thermal  resistances  of  the  Alumina  and  AIN 
substrate  and  convective  thermal  resistance  follow: 


Alumina 


0.002 

18  (0.00088)  (1) 


0.13  ’‘C/W 


AIN 


,  _  _ 0.002 

coBd  -  230  (0.00088)  (1) 


0.01  “C/R^ 
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200  (0.001)  (1)  ' 


D 

^conv 

showing  that  the  thermal  resistance  to  heat  transfer  by 
convection  is  38  and  500  times  as  large  as  resistance  to  heat 
transfer  by  conduction  in  the  Alumina  and  AlN  respectively. 

O.  EFFECT  OF  CONVECTION  HEAT  TRANSFER  ON  SYSTEM  BEHAVIOR 

This  study  focuses  on  the  effect  of  convection  cooling  on 
the  temperature  gradients  developed  within  the  low  thermal 
conductivity  substrates  and  the  overall  temperature  decrease 
of  the  chip.  The  electronic  packages  with  high  thermal 
conductivity  substrates  do  not  develop  any  temperature 
gradients  and  exhibit  an  overall  decrease  in  package 
temperature.  The  following  parameters  are  held  constant; 
volumetric  power  output  (q"*  W/m’)  of  the  chip,  ambient 
temperature  (T.)  of  the  convective  fluid,  and  the  following 
materials,  silicon  chip,  Lead/Tin  solder  and  Epoxy  Fiberglass 
substrate . 

The  governing  equation  for  convection  cooling  is  given  by: 

g"  =  h(r-r.)  (*.12) 

where  q"  is  the  heat  flux  in  W/m^  and  T  is  the  surface 
temperature.  The  inverse  relation  of  system  temperature  to 
convection  coefficient  is  shown  in  Figure  4.9.  This  figure 
illustrates  the  need  for  ever  increasing  cooling  to  maintain 
electronic  packages  within  a  temperature  tolerance  range.  An 
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CHIP  TEMPERATURE  T(*C) 


increase  in  convective  cooling  from  100  W/m^K  to  300  W/m^K 
results  in  a  reduction  of  chip  temperature  of  90  ’C.  The  same 
increase  in  convective  cooling  from  300  W/m^K  to  600  W/m^K 
produces  a  decrease  in  chip  temperature  of  only  26  *C. 


CONVECTION  COEFFICIENT  h  (W/m*«C) 


Figure  4.9  Convection  coefficient  vs.  chip  temperature 

The  temperature  gradients  in  both  X  and  Y  directions  that 
are  developed  within  the  low  thermal  conductivity  substrates 
are  diminished  by  increasing  the  convective  heat  transfer  as 
shown  in  Figures  4.10  and  4.11. 
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CONVECTION  COEFFICIENT  h  (W/m 


TEMPERATURE  GRADIENT  0,^  (®C/mm) 


Figure  4.10  9^  In  epoxy  fiberglass  substrate  vs.  convection 
coefficient 


In  the  X  direction  the  maximum  temperature  gradient  is 
decreased  by  60%  by  increasing  the  convective  cooling  by  250%. 
In  the  Y  direction  the  maximum  temperature  gradient  is 
decreased  by  52%  by  increasing  the  convective  cooling  by  250%. 


y  LOCATION  THROUGH  SUBSTRATE  (mm) 


Figure  4.11  9^  in  epoxy  fiberglass  substrate  vs.  convection 
coefficient 


E.  EFFECT  OF  GEOMETRY  OH  SYSTEM  BEHAVIOR 

The  dimensions  of  the  chip  and  substrate  are  varied  by 
assuming  the  chip  is  a  percentage  of  a  fixed  overall  height 
and  the  solder  and  substrate  make  up  the  remaining  percentage. 
We  define  the  non-dimensional  term  p  as  the  ratio  of  chip 
height  to  overall  package  height  (H)  as  in  equation  4.13, 
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(4.13) 


Substrates  are  approximately  kept  at  a  minimum  of  one  half 
a  millimeter  [Ref.  1].  The  following  material  combination  was 
selected;  a  Silicon  chip,  a  Lead/Tin  solder  and  a  Epoxy 
Fiberglass  substrate.  The  following  parameters  were  held 
constant;  convective  cooling,  ambient  temperature,  the  height 
of  the  air /solder  interface  and  the  overall  height  of  the 
package . 

First,  a  study  with  a  constant  heat  flux  (q")  was 
conducted  by  varying  the  volumetric  heat  generation  (q'")  to 
obtain  a  flux  of  4  x  10^  W/m^.  Next  the  volumetric  heat 
generation  (q'")  was  constant  which  resulted  in  different  heat 
fluxes  from  the  chip. 

1.  Constant  Heat  Flux  <q")  • 

The  result  of  maintaining  a  constant  heat  flux  is  a 
chip  that  is  a  highly  concentrated  heat  source  when  /3=0 . 18  and 
a  greatly  distributed  heat  source  when  /3=0.81.  In  Figure  4.12 
the  non-dimensional  temperature  profile  is  shown  for  0.18 
to  0.81  at  a  value  of  =  .564. 
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Figure  4.12  Hon-dlmensional  tenperature  for  various  sise  chips 
constant  heat  flux 


When  the  chip  height  results  in  a  value  of  p=0.1B,  the 
chip  is  acting  as  a  concentrated  heat  source  with  a  volvunetric 
heat  generation  q'"  =  50  MW/m^  and  the  temperature  of  the  chip 
is  greatest.  The  larger  the  chip  is  to  the  rest  of  the 
package  the  lower  the  temperature  of  the  overall  package  is. 
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A  partial  explanation  for  this  is  the  amount  of  surface  area 
exposed  to  convective  cooling  is  smaller  when  /3  =  0.18  then 
when  (3  =  0.81. 

The  thermal  strain  profile  is  shown  if  Figure  4.13  for 
the  same  value  of  (f>^  =  .564.  The  thermal  strain  profiles  can 
be  explained  in  much  the  same  way  the  temperature  profiles 
were.  That  is,  when  the  chip  is  smallest  the  largest 
expansion  of  the  package  occurs. 


Figure  4.13  Thermal  strain  profile  for  various  chip  heights 
constant  heat  flux 
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The  temperature  gradient  0^  is  shown  in  Figure  4.14 
for  /3  values  of  0.18  through  0.81.  Figure  4.14  shows  that  6 


q  =  CONSTANT 

SILICON . CHIP 

Pb  -Sn . SOLDER 

EPOXY 

FIBERGLASS-SUBSTRATE 


*  .564 


TEMPERATURE  GRADIENT  0^  (•C/mm) 


Figure  4.14  9^  for  various  chip  heights  at  constant  heat  flux 

increases  with  an  increase  in  (3,  from  =  36  *C/mm  at  /3  =  .18 
to  0y  =  57  *C/mm  at  /3  =  .81.  When  /3  =  .18  the  volumetric  heat 
generation  q'"  is  largest  at  50  MW/m^  and  provides  the  least 
surface  area  exposed  to  convection.  The  thermal  resistance 
for  the  substrate  is  e 'tributable  to  the  temperature  gradient 


5 


developed  in  the  substrate.  At  p  =  .18  the  thermal  resistance 


3.2 

.16 ( .88) 1 


«  22.7 


is  large,  thereby  not  allowing  much  heat  flux  to  conduct 
through  the  substrate,  thereby  keeping  0^  small. 

At  p  =  .81  the  thermal  resistance 

^cond(P  =  -81)  =  .16  (.^88)  1  * 

is  very  small,  allowing  for  a  much  larger  heat  flux  to  pass 
through  the  substrate  thereby  incrjasing  0^.  Figure  4.15 
shows  the  temperature  gradi'-nt  0^^  at  varying 


Figure  4.15  0,  for  various  chip  heights  at  constant  heat  flux 
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positions  of  for  (3  of  0.18  through  0.81.  These  gradients 
are  taken  at  the  top  of  the  substrate  corresponding  to  a  value 
of  =  .47.  At  =  -64  0J,  is  approximately  constant  at  97 
*C/mm  for  all  values  of  /3.  This  can  be  attributed  to  the  fact 
that  no  matter  what  thermal  energy  is  generated  by  the  chip  it 
will  tend  to  always  flow  through  the  solder  and  towards  the 
upper  right  convective  portion  of  the  substrate. 

In  the  section  of  substrate  between  chips 
corresponding  to  =  .75  and  =  .84,  increases  slightly 
from  31  *C/mm  to  52  *C/mm  and  14  *C/mm  to  27  *C/mm 
respectively.  This  can  be  accounted  for  by  looking  at  the 
substrate  equivalent  thermal  resistances.  When  /8  =  .81  the 
resistance  is  small  in  the  Y  direction  as  compared  to  the  X 
direction  as  shown  by, 


cond*  X 


1.5 

.16  (.4)  1 


23.4®C/W 


^  ^cond> 


y 


.4 

.16  (.5)  1 


5°C/W 


Therefore  heat  transfer  would  rather  take  place  in  the  Y 
direction  instead  of  the  X  direction.  When  /3  =  .18,  the 
thermal  resistances  are  large  in  the  Y  direction  and  small  in 
the  X  direction  that  is. 
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■  .16  (3,2)  1  “  3°C/W 

'  7iT(.Vt  ' 

This  change  in  thermal  resistance  now  favors  heat  transfer  in 
the  X  direction,  which  accounts  for  the  increase  in  0,,  for 
smaller  values  of  /3. 

2.  Constant  Volumetric  Heat  Generation  (g'")  • 

This  study  shows  the  effect  of  using  larger  power 
chips  in  increasing  sizes.  The  heat  generation  of  each  chip 
is  proportional  to  chip  volume,  that  is,  heat  generation 
increases  with  an  increase  in  /3.  Figure  4.16  is  the 


Figure  4 . 16  Mon-dlRenslonal  temperature  for  various  size  chips 
constant  volumetric  power  generation 
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normalized  temperature  profile  for  (3  =  0.18  to  ft  =  0.81  at  a 
position  of  =  .564. 

The  non  dimensional  temperature  of  the  package  follows 
the  identical  relationship  as  did  volumetric  heat  generation 
to  p.  As  /3  increases  so  does  the  non-dimensional  temperature 
of  the  package.  If  the  ratio  of  heat  flux  to  non-dimensional 
temperature  is  examined,  the  increase  in  non-dimensional 
temperature  is  found  to  be  non-linear.  The  following  ratios 
q'7T^  are  given  for  various  values  of  p, 


p 

.  18 

m 

.60 

.73 

.81 

j/f 

f  w/mM 

10 

11.2 

16.1 

16.0 

16.0 

The  increase  in  heat  flux  q"(W/m^)  from  p  =  .18  to  ^  =  .31 
is  12  KW/m^  and  from  /3  =  .6  to  /3  =  .81  the  increase  is  14 
KW/m^.  Even  though  the  increase  in  q"  is  2  KW/m^  more  for  the 
later  /3's  the  ratio  of  q'VT^^  remains  constant.  The  factors 
that  contribute  to  this  are,  an  increase  in  the  convection 
surface  area  of  the  chip  and  a  decrease  in  the  conduction 
thermal  resistance  of  the  substrate  to  allow  for  a  greater 
heat  transfer  through  the  substrate. 

Figure  4.17  shows  the  thermal  strain  profile  at  a  value  of 
=  .564.  The  chip  does  not  have  a  large  increase  in  thermal 
strain  over  the  range  of  /3's  while  the  solder  has  a  very  large 
increase  in  thermal  strain.  This  goes  back  to  Table  4.1, 
where  the  coefficient  of  thermal  expansion  for  the  Pb-Sn 
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solder  is  roughly  eleven  times  greater  than  that  of  Silicon. 


THERMAL  STRAIN  €.=  ( T -T„  )  a  (  X  lo’  m/m) 

T  ® 


Figure  4.17  Thermal  strain  profile  for  various  chip  heights 
constant  volumetric  power 

Figure  4.18  shows  the  variation  of  the  temperature 
gradient  9^  at  the  top  of  the  substrate  at  various  values  of 
At  values  of  =  0.64  and  -  0.43  the  temperature 
gradient  is  largest  and  this  corresponds  to  either  side  of  the 
solder  joint.  At  other  values  of  the  temperature  gradient 
is  diminished  greatly.  At  a  value  of  ft  =  0.81  the  temperature 
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gradients  are  large.  They  are  almost  twice  as  large  as 


encountered  in  any  other  analysis.  For  this  large  a  chip 
(large  q")  the  convection  cooling  must  be  increased 
significantly  to  reduce  this  temperature  gradient.  For  small 
chip  sizes  the  temperature  gradient  is  very  small  for  any 
value  of  <t>^. 


Figure  4.18  6,  for  various  chip  heights  at  constant  voluamtric 
power 


5 


Figure  4.19  shows  the  variation  in  0^  at  the  top  of 
the  substrate  where  the  largest  value  of  the  temperature 
gradient  within  the  depth  of  the  substrate  occurs . 

When  p  =  0.81  the  chip  has  the  greatest  heat  flux, 
q"=72  KW/m^  with  a  correspondingly  large  0^.  When  )3  =  0.18  the 
chip  has  the  lowest  heat  flux,  q'^  16  KW/m^  and  0^  is  all  but 
non-existent . 


TEMPERATURE  GRADIENT  ©  (•C/mm) 

Figure  4.19  6^  for  various  chip  heights  at  constant  voluaetric 
l^wer 
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3.  Overview  Of  OeMietric  Study 

For  both  the  constant  heat  flux  q"  (W/m^)  and  the 
constant  volumetric  heat  generation  q'"  (W/m^),  the  temperature 
gradient  in  the  Y  direction  at  the  top  of  the  substrate  is 
greatest  when  the  chip  is  largest  (jS  =  .81).  The  reasons  as 
given  previously  are  the  interaction  between  the  conductive 
thermal  resistance,  the  convective  thermal  resistance,  and  the 
ratio  of  volumetric  heat  generation  to  available  convection 
area. 

If  the  convective  thermal  resistance  is  small  as 
compared  to  the  conductive  thermal  resistance  then  the 
convective  mode  will  dominate  the  heat  transfer  process  as  the 
heat  transfer  process  is  governed  by  the  path  of  least 
resistance. 

F.  EFFECT  OF  CHIP  ENCAPSUIATIOH  ON  SYSTEM  TEMPERATURE 

Some  electronic  packages  are  protected  from  moisture  and 
the  environment  by  a  protective  "encapsulating"  coating.  For 
this  study  the  encapsulation  material  EME-llOO-T. [Ref .  1]  is 
0.1875  mm  thick.  The  effect  of  encapsulation  on  the  chip 
temperature  and  gradients  developed  within  the  package  is 
studied. 

The  following  parameters  will  be  held  constant,  volumetric 
heat  generation  (q'")#  convection  cooling  (h),  and  ambient 
temperature  (T,).  These  values  were  the  same  as  in  section  B 
for  comparison  purposes.  The  materials  are  a  Silicon  chip,  a 
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Lead/Tin  Solder,  with  all  four  substrates,  and  the 
encapsulation  material.  The  thermal  conductivity  of  the 
encapsulation  material  is  k  =  0.67  W/m*C.  Within  the  other 
encapsulation  materials  the  thermal  conductivity  (k)  varied 
from  0.67  to  1.97  W/m’C. 

Figure  4.20  is  the  non-dimensional  temperature  profile 
taken  at  a  value  of  =  0.564.  The  entire  package 
experiences  less  then  a  5%  temperature  change,  when  low 
thermal  conductivity  substrates  are  used  and  less  then  2% 
temperature  change  when  high  thermal  conductivity  substrates 
are  used.  The  encapsulation  material  causes  the  surface 
temperature  to  decrease  by  a  small  amount. 


T-T„ 

NON  DIMENSIONAL  TEMPERATURE  T  » 


Figure  4.20  Encapsulated  chip  non-diaenaional  teaperature 


NON  DIMENSIONAL  Y  LOCATION  <f> 


Even  though  the  encapsulation  material  has  a  low  thermal 
conductivity  its  overall  thermal  resistance  is  low  because  its 


less  then  two  tenths  of  a  millimeter  thick. 

R  .  =  - ■  187  5  11°C/W 

.67(2.5)1 

This  is  a  very  low  thermal  resistance  accounting  for  the  small 
change  in  temperature  profile.  Figure  4.21  shows  the  thermal 
strain  profile  for  <f>^  =  .564. 


Alumina  epoxy  fiberolass 


Figure  4.21  Encapsulated  chip  thermal  strain  profile 
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The  thermal  strain  profile  has  the  same  negligible 
increases  as  the  normalized  temperature  profile  had  except  the 
encapsulation  material  has  a  very  large  coefficient  of  thermal 
expansion  (a). 

In  summary  the  effect  of  encapsulation  on  system 
temperature  and  temperature  gradients  are  negligible. 
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V.  STRESS  AMALXSIS 


A.  OVERVIEW  OP  STRESS  ANALYSIS 

The  stress  analysis  conducted  is  based  on  using  different 
materials  in  the  package  and  keeping  all  other  factors 
constant.  A  silicon  chip  material  was  used  for  all  cases, 
with  two  solders  Pb-Sn  and  Au-Sn  ,  and  the  four  substrates 
Epoxy  Fiberglass,  Polyimide  Fiberglass,  Alumina  and  Aluminum 
Nitride.  These  eight  sets  of  material  combinations  correspond 
to  the  thermal  strain  graphs,  shown  in  Figures  4.3  and  4.4  in 
Chapter  IV. 

Figure  4 . 3  corresponds  to  the  material  combinations  of  a 
silicon  chip,  the  Pb-Sn  solder,  and  the  four  substrates  listed 
above.  Figure  4.4  corresponds  to  the  material  combinations  of 
a  silicon  chip,  the  Au-Sn  solder,  and  the  four  substrates 
listed  above. 

The  following  other  factors  remain  constant,  q"'=  20  MW/m^ 
and  h  =  250  W/m^*C.  All  stresses  given  with  the  exception  of 
the  bending  stresses  at  the  centerline  oc  the  package  are 
along  the  solder  interfaces.  The  following  stresses  are 
calculated.  First,  the  bending  stresses  (a^,),  along  the 
package  centerline.  Three  stresses  are  calculated  along  the 
solder-chip  interface,  the  upper  shear  stress  (r^),  an  upper 
bending  stress  (cr^^u)/  and  an  upper  normal  stress  (o^  u).  Three 
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stresses  are  calculated  along  the  solder-substrate  interface, 
the  lower  shear  stress  (t^),  a  lower  bending  stress  l)/  and 
a  lower  normal  stress 

Table  5 . 1  gives  mechanical  properties  of  electronic 
packaging  materials,  elastic  modulus  given  in  GPa.  These 
properties  are  used  in  the  calculation  of  stresses  in 
accordance  with  the  FEM  development  of  chapter  three. 

TABLE  5.1 


Material 

Elastic  Modulus 

Poisson's  Ratio 

Silicon 

13.03 

0.3 

Pb-Sn 

7.4 

0.4 

Au-Sn 

59.2 

0.4 

EP  Fbgls 

35.0 

0.3 

Polyimide 

40.0 

0.3 

Alumina 

262.0 

0.3 

AIN 

339.0 

0.3 

B.  STRESS  AHALTSIS  WITH  (Pb-Sn)  SOLDER 
1.  Bending  Stress  Along  Centerline, 

The  bending  stress  at  the  centerline  of  the  package  is 
the  lowest  of  all  stresses  presented  and  is  shown  in  Figure 
5.1.  The  bending  stress  a^,  through  the  chip  is  of  the  order 
of  1  MPa  for  all  four  substrates.  Through  the  chip  the 
bending  stress  is  tensile  for  the  low  thermal  conductivity 
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substrates  (Epoxy  and  Polyimide  Fiberglasses)  and  compressive 
for  the  high  thermal  conductivity  substrates  (Alumina  and 
AlN) .  Through  the  substrates  the  bending  stresses  in  the  low 


-10  -8  -6  -4  -2  0  2  4  6 


BENDING  STRESS  (t  MPc) 

Figure  5 . 1  Centerline  bending  stresses  a,, 


thermal  conductivity  substrates  is  a  maximum  in  compression  at 
the  top  of  the  substrate  =  .47)  at  -8  MPa  and  is  a  maximum 
in  tension  at  the  bottom  of  the  substrate  (0^  =  0)  at  6  MPa. 
The  high  thermal  conductivity  substrates  have  just  the 
opposite  results,  with  a  maximum  in  tension  at  the  top  of  the 
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substrate  {<py  =  .47)  at  2.5  Mpa  and  maximum  in  compression  at 
the  bottom  of  the  substrate  ((p^  -  0)  at  -2.5  MPa.  Thus  the 
bending  stresses  in  the  low  thermal  conductivity  substrates 
are  about  twice  as  large  as  the  bending  stresses  in  the  high 
thermal  conductivity  substrates.  The  advantage  with  respect 
to  bending  stresses  along  the  centerline  of  the  chip  is  with 
the  high  thermal  conductivity  substrates  (alumina  and  AlN). 

2.  Upper  Nonaal  Stress, 

The  upper  normal  stress  distribution  along  the  chip- 
solder  interface  for  all  the  substrates  are  compressive  with 
a  minimum  occurring  at  6  =  .05  and  maximum  occurring  at  6  ® 
0.8  as  shown  in  Figure  5.2.  For  the  low  thermal  conductivity 
substrates  the  maximum  upper  normal  stress  for  Epoxy 
Fiberglass  is  -49  MPa,  and  the  maximum  for  Polyimide 
Fiberglass  is  -40  Mpa.  Both  curves  show  a  decrease  in  the 
upper  normal  stress  at  the  edges  of  the  solder  due  to  the  fact 
that  the  edges  at  6  =  0.0  and  6=1.0  are  free  surfaces.  For 
the  Pb-Sn  group  of  materials  the  normal  stress  is  the  maximum 
stress  examined  for  the  low  thermal  conductivity  substrates. 

The  upper  normal  stresses  for  the  high  thermal 
conductivity  substrates  are  about  -16  MPa  and  are  within  5%  of 
each  other,  decreasing  in  magnitude  slightly  at  the  edges  of 
the  solder.  This  shows  that  the  upper  normal  stresses  in  the 
high  thermal  conductivity  substrates  are  about  a  third  of  the 
upper  normal  stress  in  the  low  thermal  conductivity 
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substrates.  We  note  that  the  normal  stresses  are  fairly 
uniform  from  6  -  .25  to  6  =  .9. 

In  contrast  to  the  chip-solder  interface,  the  normal 
stresses  at  the  solder-substrate  interface  (i.e.  the  lower 
normal  stresses)  are  very  small  on  the  order  of  KPa  and  in 
addition  are  tensile  stresses  for  all  substrate  combinations. 
With  respect  to  ^  the  advantage  is  again  seen  with  the  high 
thermal  conductivity  substrates  (alumina  and  Ain). 


0  .1  .2  .3  .4  .5  .6  .7  .8  .9  1.0 

NORMALIZED  X  LOCATION  ALONG  SOLDER  8  ■ 

L 


Figure  5.2  Upper  normal  stresses 
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3.  Upper  Bending  Stress/ 

The  upper  bending  stress  ^  along  the  chip  solder 
interface,  is  defined  as  the  normal  stress.  However,  these 
stresses  are  presented  along  the  interface,  where  they  are 
greatest,  as  a  function  of  6  rather  than  across  the  Y 
direction  through  the  solder.  These  bending  stresses  are 
compressive  for  all  substrate  combinations.  For  the  low 
thermal  conductivity  substrates  the  maximum  upper  bending 
stresses  occur  at  6  =  0.05  and  0.95  with  values  of  -16  MPa  and 
-15  Mpa  respectively.  For  the  high  thermal  conductivity 
substrates  maximum  upper  bending  stresses  occur  at  6  =  0.05 
with  -17  Mpa  magnitude,  and  minimums  of  -8.5  MPa  occur  at  6  = 
0.8  as  shown  in  Figure  5.3.  Again  we  note  that  severe 
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gradients  occur  close  to  the  free  edges,  and  with  fairly 
uniform  stresses  along  the  center  of  the  chip-solder 
interface.  We  note  these  bending  stresses  fall  between  the 
bending  stresses  at  the  centerline  and  normal  stresses 
reported  in  the  previous  sections.  With  respect  to  these 
stresses,  the  advantage  (i.e.,  lower  stresses)  is  obtained  for 
the  low  thermal  conductivity  substrates  (epoxy  and  polyimide 
fiberglasses). 

4.  Lower  Bending  Stress 

The  lower  bending  stresses  are  compressive  in  nature 
for  all  substrate  combinations  as  shown  in  Figure  5.4.  For 
the  low  thermal  conductivity  substrates,  an  immediate  decrease 
from  about  -18  MPa  to  about  -13  MPa,  that  is  about  a  30% 
change.  Thereafter  a  gradual  decrease  to  about  -11  MPa,  that 
is,  about  15%  change  over  the  remaining  domain.  The  lower 
bending  stress  for  the  high  thermal  conductivity  substrates 
remains  constant  over  the  entire  length  of  the  solder  joint  at 
-15  MPa. 

In  comparison  with  the  upper  bending  stresses  the 
lower  bending  stresses  are  greater  for  the  low  thermal 
conductivity  substrates.  The  advantage  for  the  lower  bending 
stresses  goes  with  the  high  thermal  conductivity  substrates 
(alumina  and  AlN) . 
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LOWER  BENDING  STRESS  <r  (MPa) 
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Figure  5.4  Lower  bending  stress  0|,  | 

5.  Upper  Shear  Stress 

The  upper  shear  stresses  at  the  chip/solder  interface 
are  positive  at  6  =  .05  and  negative  at  6  *  .95  for  all 
substrate  combinations  as  shown  in  Figure  5.5.  For  the  low 
thermal  conductivity  substrates  the  maximum  shear  stress  is 
negative  and  occurs  at  6  =  .95  at  -14  MPa  and  crosses  zero 
stress  at  6  »  .15.  For  the  high  thermal  conductivity 
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substrates  the  maximum  shear  stress  is  positive  and  occurs  at 
6  =  .05  at  11  MPa  and  crosses  zero  stress  at  6  =  .8.  The 
upper  shear  stresses  have  large  gradients  at  either  edge  of 
the  chip-solder  interface,  but  remain  fairly  constant  through 
the  majority  of  the  domain.  From  6  =  .2  through  6  =  .9  the 
upper  shear  stresses  are  less  than  5  MPa  in  magnitude. 

These  stresses  fall  in  between  the  upper  bending 
stresses  and  the  bending  stresses  at  the  centerline.  With 
respect  to  these  stresses  the  advantage  goes  to  the  high 
thermal  conductivity  substrates  (alumina  and  AlN) . 
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Figure  5.5  Upper  shear  stress 
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6.  Lower  Shear  Stress 

The  lower  shear  stresses  are  opposite  in  sign  from  the 
upper  shear  stresses,  in  that  the  lower  shear  stresses  at  the 
solder/substrate  interface  are  negative  at  6  =  .05  and  are 
positive  at  5  =  .95  for  all  substrate  combinations  as  shown 
in  Figure  5.6.  For  the  low  thermal  conductivity  substrates 
the  maximum  shear  stresses  occur  at  6  =  .05  at  about  -16  MPa. 
This  is  comparable  in  magnitude  with  the  upper  shear  stresses. 
For  the  high  thermal  conductivity  substrates  the  maximum 
positive  shear  stress  occurs  at  6  =  .95  at  16  MPa  while  the 
maximum  negative  shear  stress  occurs  at  6  =  .05  at  -10  Mpa. 
This  stress  at  6  =  .05  is  twice  the  magnitude  as  compared  to 
the  upper  shear  stresses.  Again  we  note  that  steep  gradients 
occur  at  the  edges  of  the  solder-substrate  interface  and  the 
low  thermal  conductivity  substrates  have  a  gradual  decrease 
over  the  center  of  the  domain.  The  stresses  of  the  low 
thermal  conductivity  substrates  cross  zero  at  6  =  .9.  The 
stresses  for  the  high  thermal  conductivity  substrates  vary 
linearly  with  6  between  6  =  .2  and  .9.  These  stresses 
range  from  -5  MPa  at  6  =  .15  to  10  Mpa  at  6  =  .9  crossing  zero 
stress  at  6  =  .45.  We  note  these  lower  shear  stresses  are 
comparable  in  magnitude  to  the  upper  shear  stresses.  With 
respect  to  these  stresses,  there  is  no  clear  advantage  for 
either  set  of  substrates. 
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Figure  5.6  Lower  shear  stress 


C.  STRESS  ANALYSIS  WITH  (Au-Sn)  SOLDER 

1.  Bending  Stress  Along  Centerline,  o,, 

The  bending  stress  at  the  centerline  of  the  package 
remains  the  lowest  of  all  stresses  presented  as  shown  in 
Figure  5.7.  The  bending  stresses  at  the  centerline  =0.0) 
for  the  packages  with  Au-Sn  solder  are  identical  to  those  with 
Pb-Sn  solder.  The  centerline  bending  stress  is  independent  of 
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which  solder  is  used.  With  respect  to  these  stresses,  the 
advantage  is  obtained  by  the  high  thermal  conductivity 
substrates  (alumina  and  AlN). 


Figure  5.7  Centerline  bending  stress  o,, 

2.  Upper  Normal  Stress, 

The  upper  normal  stress  distribution  along  the  chip- 
solder  interface  for  all  the  substrates  are  compressive  in 
nature  as  shown  in  Figure  5.8.  For  the  low  thermal 
conductivity  substrates  the  maximum  for  Epoxy  Fiberglass  is 
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-46  MPa  and  the  maximum  for  Polyimide  Fiberglass  is  -40  Mpa 
and  they  are  constant  over  a  range  of  .  1  <  6  <  .  9 .  Both  upper 


Figure  5.8  Upper  normal  stress  o, 


o,u 


normal  shear  stress  distributions  decrease  at  the  edges  of  the 
solder  due  to  the  fact  that  the  stress  at  a  free  surface  is 
zero.  The  normal  stress  is  once  again  the  maximum  stress 
examined  for  the  low  thermal  conductivity  substrates. 

The  high  thermal  conductivity  substrates  have 
increased  by  about  20%  over  the  Pb-Sn  material  package.  There 
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are  sharp  increases  at  5  -  .1  and  a  slight  increase  to  a 
maximum  of  about  -20  MPa  at  6  =  .9. 

Again  we  note  that  in  contrast  to  the  chip-solder 
interface,  the  normal  stresses  at  the  solder  substrate 
interface  are  very  small  on  the  order  of  KPa  and  in  addition 
are  tensile  stresses  for  all  substrate  combinations.  With 
respect  to  these  stresses,  the  advantage  is  obtained  for  the 
high  thermal  conductivity  substrates  (alumina  and  AIM)  and 
with  respect  to  solders,  the  advantage  is  obtained  for  the  Pb- 
Sn  solder. 

3.  Upper  Bending  Stress,  Oj,,, 

The  upper  bending  stress  ^  along  the  chip-solder 
interface  is  defined  as  the  normal  o^  stress  as  previously 
mentioned.  The  upper  bending  stresses  are  compressive  for  all 
substrate  combinations  as  shown  in  Figure  5.9.  For  the  low 
thermal  conductivity  substrates  maximums  occur  at  6  =  0.05  and 
0.95  with  values  of  -19  MPa  and  -22  Mpa  respectively.  For  the 
high  thermal  conductivity  substrates  a  maximum  occurs  at  6  = 
0.05  at  -39  Mpa.  This  upper  bending  stress  is  twice  as  large 
for  the  high  thermal  conductivity  substrates  with  Au-Sn  solder 
then  Pb-Sn  solder.  There  is  a  very  large  gradient  at  6  =  0.05 
where  the  upper  bending  stress  decreases  from  -39  MPa  to  -22 
MPa,  that  is,  about  a  45%  decrease.  For  the  high  thermal 
conductivity  substrates  the  upper  bending  stress  continues  to 
decrease  at  a  gradual  rate  across  the  remaining  domain,  6  =  .2 
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to  6  =  .9.  With  respect  to  these  stresses,  the  advantage  is 
obtained  for  the  low  thermal  conductivity  substrates  (epoxy 
and  polyimide  fiberglasses)  and  with  respect  to  solders,  the 
advantage  is  obtained  for  the  Pb-Sn  solder. 


Figure  5.9  Upper  Bending  Stress  ai,^ 

4.  Lower  Bending  Stress,  Oj,  ^ 

The  lower  bending  stresses  are  compressive  for  all 
substrate  combinations  as  shown  in  Figure  5.10.  For  the  low 
thermal  conductivity  substrates  a  maximum  occurs  at  6  =  0.05 


with  a  value  of  -39  MPa  and  a  minimum  at  -7  MPa  at  6  =  0.9. 


This  is  an  increase  of  200%  over  the  same  package  with  Pb-Sn 
solder.  The  lower  bending  stress  for  the  high  thermal 
conductivity  substrates  remain  fairly  constant  over  a  range  of 
.15  <  6  <  .8  with  a  value  of  -38  MPa  and  is  a  maximum  at  6  = 
.95  at  -49  MPa.  This  is  an  increase  of  about  260%  over  the 
same  package  with  Pb-Sn  solder.  The  lower  bending  stresses 
are  the  highest  stresses  presented  for  the  high  thermal 
conductivity  substrates.  With  respect  to  these  stresses,  the 
advantage  is  obtained  for  the  low  thermal  conductivity 
substrates  (epoxy  and  polyimide  fiberglasses)  and  with  respect 
to  solders,  the  advantage  is  obtained  for  the  Pb-Sn  solder. 


Figure  5 . 10  Lower  bending  stress  0|,  ^ 


5.  Upper  Shear  Stress 

The  upper  shear  stresses  at  the  chip/solder  interface 
are  positive  at  6  =  .05  and  negative  at  6  =  .95  for  all 
substrate  combinations  as  shown  in  Figure  5.11.  For  the  low 


Figure  5.11  Upper  shear  stress 

thermal  conductivity  substrates  the  maximum  shear  stress  is 
negative  and  occurs  at  6  =  .95  with  a  value  of  -15  MPa  and 
crosses  zero  stress  at  6  »  .15.  For  the  high  thermal 
conductivity  substrates  the  maximum  shear  stress  is  positive 
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and  occurs  at  6  =  .05  with  a  value  of  16  MPa  and  crosses  zero 
stress  at  6  »  .8.  The  maximum  upper  shear  stresses  for  both 
types  of  substrates  show  only  a  slight  increase  over  the  Pb-Sn 
solder  packages.  With  respect  to  these  stresses,  there  is  no 
clear  cut  advantage  for  either  set  of  substrates,  nor  either 
solder. 

6.  Lower  Shear  Stress 

The  lower  shear  stresses  are  opposite  in  sign  from  the 
upper  shear  stresses,  in  that  the  lower  shear  stresses  at  the 
solder/substrate  interface  are  negative  at  6  =  .05  and  are 
positive  at  5  =  .95  for  all  substrate  combinations  as  shown  in 
Figure  5.12.  For  the  low  thermal  conductivity  substrates  the 
maximum  shear  stresses  occur  at  6  =  .05  with  a  value  of  -21 
MPa.  This  represents  an  increase  of  130%  over  the  same 
package  with  Pb-Sn  solder.  The  lower  shear  stress  then 
decreases  to  -10  Mpa  at  6  =  .1  and  then  slowly  decreases  to 
about  zero  at  5  =  .95.  For  the  high  thermal  conductivity 
substrates  the  maximum  shear  stresses  occur  at  6  =  .95  with  a 
value  of  39  MPa.  This  represents  an  increase  of  240%  over  the 
same  package  with  Pb-Sn  solder.  At  6  =  .05  the  lower  shear 
stress  is  negative  at  -32  MPa  which  represents  an  increase  of 
110%  over  the  ssune  package  with  Pb-Sn  solder.  With  respect  to 
these  stresses,  the  advantage  is  obtained  for  the  low  thermal 
conductivity  substrates  (epoxy  and  polyimide  fiberglasses)  and 
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LOWER  SHEAR  STRESS  r.  (MPa) 


with  respect  to  solders,  the  advantage  is  obtained  for  the  Pb 
Sn  solder. 
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Figure  5.12  Lower  Shear  Stress 

D.  SUMMARY  OF  STRESS  ANALYSIS 

The  following  table  represents  the  maximum  stress  '  (given 
in  MPa)  for  each  material  combination  studied.  Various 
comparisons  are  possible. 
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EP  FBGLS 


POLY  FBGLS 


ALUMINA 


AIN 


TABLE  5.2 

MAXIMUM  STRESSES 

Pb-Sn 

Au-Sn 

-8 

-8 

-49 

-46 

-16 

-20 

<^b,l 

-19.5 

-39 

-14.5 

-16 

-16 

-21 

‘^b 

-7.5 

-7.5 

-41 

-40 

*^b,d 

-14 

-18 

^b,l 

-16.5 

-34 

-13 

-14 

-15 

-20 

^b 

-2.5 

-3 

-16 

-20 

^b,u 

-17 

-38 

^b,l 

-15 

-45 

11 

16 

15 

40 

‘^b 

-2.5 

-3 

-15 

-19 

‘^b,u 

-17 

-38 

^b,l 

-14 

-43 

11 

16 

15 

40 
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For  the  Epoxy  fiberglass  substrate  packages  there  is 
little  difference  in  centerline  bending  a^,  normal  and 
upper  shear  stresses  and  only  a  slightly  larger  upper 
bending  and  lower  shear  Tj  stresses  for  Au-Sn  versus  Pb- 
Sn  solder.  For  this  substrate  the  lower  bending  ^  stresses 
are  twice  as  great  for  Au-Sn  versus  Pb-Sn  solder.  The  normal 
stresses  are  the  largest  stresses  for  this  substrate  for  both 
Pb-Sn  and  Au-Sn  solders.  The  most  significant  stresses  with 
regard  to  deforming  the  solder  joint  are  the  lower  bending  and 
lower  shear.  The  upper  normal  stress  is  not  significant 
because  it  tends  to  compress  a  solder  crack  rather  then 
propagate  it. 

For  the  Polyimide  fiberglass  substrate  packages  there  is 
little  difference  j.t-.  -enterline  bending  normal  and 
upper  shear  stresses  and  only  a  slightly  larger  upper 
bending  Oj,  and  lower  shear  stresses  for  Au-Sn  versus  Pb- 
Sn  solder.  For  this  substrate  the  lower  bending  ^  stresses 
are  twice  as  great  for  Au-Sn  versus  Pb-Sn  solder.  The  normal 
stresses  are  the  largest  stresses  for  this  substrate  for  both 
Pb-Sn  and  Au-Sn  solders.  The  most  significant  stresses  with 
regard  to  deforming  the  solder  joint  are  the  lower  bending  and 
lower  shear.  The  upper  normal  stress  is  not  significant 
because  it  tends  to  compress  a  solder  crack  rather  then 
propogate  it.  The  Polyimide  fiberglass  packages  have  stress 
distributions  similiar  to  Epoxy  fiberglass,  but  tend  to  be 
slightly  lower  in  magnitude. 
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For  the  Alumina  substrate  only  the  centerline  bending 
stresses  are  similiar  for  both  solders.  The  upper  normal 
and  upper  shear  stresses  are  slightly  larger  for  Au-Sn 
versus  Pb-Sn  solders.  The  upper  bending  ^  stress  is  225% 
larger,  lower  bending  stress  is  300%  larger,  and  the  lower 
shear  stress  is  265%  larger  for  the  Au-Sn  versus  Pb-Sn 

solder.  The  upper  bending  stress  is  largest  for  the  Pb-Sn 

solder  while  the  lower  bending  stress  is  largest  for  the  Au-Sn 
solder.  The  most  significant  stresses  with  regard  to 
deforming  the  solder  joint  are  the  upper  bending,  lower 
bending  and  lower  shear. 

For  the  AIN  substrate  only  the  centerline  bending  a^^ 
stresses  are  similiar  for  both  solders.  The  upper  normal 
and  upper  shear  stresses  are  slightly  larger  for  Au-Sn 
versus  Pb-Sn  solders.  The  upper  bending  stress  is  225% 
larger,  lower  bending  Oj,  ^  stress  is  305%  larger,  and  the  lower 
shear  stress  is  265%  larger  for  the  Au-Sn  versus  Pb-Sn 

solder.  The  upper  bending  stress  is  largest  for  the  Pb-Sn 

solder  while  the  lower  bending  stress  is  largest  for  the  Au-Sn 
solder.  The  most  significant  stresses  with  regard  to 
deforming  the  solder  joint  are  the  upper  bending,  lower 
bending  and  lower  shear.  The  AlN  packages  have  stress 
distributions  which  are  approximately  less  than  5%  different 
from  Alumina. 

Based  on  the  stress  results  from  this  investigation  the 
following  ranking  of  material  combinations  are  presented  based 
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on  obtaining  the  lowest  stresses  in  the  package.  The  silicon 
chip  is  used  in  all  packages. 


1.  Pb-Sn  with  Altunina  or  Ain 

2.  Pb-Sn  with  Polyimide  fiberglass 

3.  Pb-Sn  with  Epoxy  fiberglass 

4.  Au-Sn  with  Polyimide  fiberglass 

5.  Au-Sn  with  Epoxy  fiberglass 

6.  Au-Sn  with  Alumina  or  Ain 

Bending  (aj,)  and  Shear  (t)  stresses  at  the  solder- 
substrate  interface  can  be  significant.  The  lower  inside 
corner  of  the  solder-substrate  interface  (5  =  0.0)  is  the 
point  where  failure  is  most  likely  to  occur  due  to  bending  and 
shear  stresses. 

This  investigator  recommends  the  following  topics  for 
further  research.  Alter  the  thermal  code  for  transient 
temperature  response  of  the  package.  The  result  of  this  would 
be  to  track  stress  in  time.  In  this  way  the  changes  in 
temperature  and  stress  can  be  tracked  in  time  after  the 
electronic  device  has  been  turned  on.  Conduct  a  parametric 
study  of  a  general  tri-material  package  varying  the  materials 
coefficient  of  thermal  expansion  and  modulus  of  elasticity  to 
better  understand  the  effect  on  stresses  produced  in  the  tri¬ 
material  package. 
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APPENDIX  A. 


PROGRAM  HEAT 

********************************************************************** 
*  * 

*  This  F.E.M.  program  solves  for  the  Temperature  distribution  in  * 

*  electronic  packages.  The  governing  differential  equation  is  * 

*  POISSON'S  equation.  More  specifically:  * 

*  * 

*  * 

*  When  symmetry  is  envoked  homogeneous  Neumann  boundary  conditions  * 

*  are  used  at  these  locations.  The  other  boundaries  are  Cauchy.  * 

*  * 

*  * 

*  Listing  of  variables:  * 

*  * 

*  aaij  -  local  nodal  point  values  of  BIGA  matrix  * 

*  ALPHA 

*  AREA  -  vector  of  element  areas  * 

*  bi  -  triangular  coordinate  property  of  each  element  * 

*  bbij  -  contribution  from  Cauchy  boundary  to  BIGA  matrix  * 

*  BIGA  -  global  left  hand  matrix  * 

*  BIGF  -  right  hand  side  vector  * 

*  BDRYVAL  -  value  of  boundary  point  * 

*  BLOW  -  (BUP,BVERT)  inverse  of  the  absolute  temp  of  the  fluid  * 

*  cauchfi  -  contribution  from  Cauchy  boundary  to  BIGF  vector  * 

*  ci  -  triangular  coordinate  property  of  each  element  * 

*  DELI, 2, 3-  difference  between  assumed  and  calc,  surface  temps  * 

*  DELMAX  «  maximum  of  dell  or  del2  or  del3  * 

*  DIVSUB  *  #  of  rows  of  substrate  * 

*  DIVX  -  #  of  divisions  in  x  direction  for  a  grid  density  * 

*  DIVY  -  #  of  divisions  in  y  direction  for  a  grid  density  * 

*  DIVXAIR  -  #  of  columns  of  air  in  grid  * 

*  DIVYAIR  -  #  of  rows  of  air( solder)  in  grid  * 

*  ffl  -  contribution  of  source  term  to  BIGF  vector  * 

*  GEN  ■  if  equal  to  1  then  internal  heat  generation  is  involved  * 

*  HLOW  -  lower  side  heat  transfer  coefficient  * 

*  HUP  -  upper  side  heat  transfer  coefficient  * 

*  HVERT  -  right  side  heat  transfer  coefficient  * 

*  IMAT  ■  identifies  what  material  properties  go  with  which  element* 

*  ICOR  -  establishes  correspondence  between  local  and  global  coor.* 

*  KEL  ■  therm,  conductivity  of  element  * 

*  KlNVS  -  kinematic  viscosity  of  the  fluid  * 

*  NBCLOW  -  if  «  1  ;  then  lower  boundary  is  a  Cauchy  * 

*  NBCUP  ■  if  -  1  ;  then  upper  boundary  is  a  Cauchy  * 

*  NBCVERT  -  if  -  1  ;  then  right  boundary  is  a  Cauchy  * 

*  NBDLEFT  -  if  «  1  ;  then  left  boundary  is  a  Direchlet  * 

*  NBDRIT  -  if  -  1  ;  then  right  boundary  is  a  Direchlet  * 

*  NBDUP  -  if  -  1  ;  then  upper  boundary  is  a  Direchlet  * 

*  NBDLOW  -  if  -  1  ;  then  lower  boundary  is  a  Direchlet  * 

*  NCOL  -  number  of  columns  for  grid  subroutine  * 

*  NROW  -  number  of  rows  for  grid  subroutine  * 

*  NSUBEL  -  #  of  substrate  elements  * 

*  NAIREL  -•  #  of  columns  of  air  in  grid  * 

eo 


*  NARSOEL 

*  QDOT 

*  HALLOW 

*  RALUP 

*  RALVERT 

*  RHO 

*  RHOEL 

*  SPHT 

*  SPHTEL 

*  TAMB 

*  TEMP 

*  TFLLOW 

*  TFLUP 

*  TFLVERT 

*  TSLOW 

*  TSUP 

*  TSVERT 

*  THERCON 

*  X 

*  XPOS 

*  Y 

*  YPOS 


-#  of  air  &  solder  elements 

-  volumetric  heat  generation  of  electronic  device 

-  Raleigh  #  of  the  fluid  at  the  lower  surface 

-  Raleigh  #  of  the  fluid  at  the  upper  surface 

-  Raleigh  #  of  the  fluid  at  the  vertical  surface 

-  density  of  4  different  materials 

-  density  of  element 

-  specific  heat  of  4  different  materials 

-  specific  heat  of  element 
ambient  temperature 

“  solution  vector  of  temperatures 

-  fluid  temperature  at  lower  surface 

■  fluid  temperature  at  upper  surface 

-  fluid  temperature  at  vertical  surface 

-  average  surface  temperature  of  lower  surface 

-  average  surface  temperature  of  upper  surface 

■  average  surface  temperature  of  vertical  surface 

■  thermal  conductivity  of  4  different  materials 

■  X  position  of  system  nodal  point 

-  X  position  of  new  grid  density 

-  y  position  of  system  nodal  point 
“  y  position  of  new  grid  density 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


INCLUDE  'HEATCOMN.FOR' 

********  ensure  the  proper  amount  of  workspace  is  allocated  *** 
********  on  problems  with  a  large  quantity  of  unknowns  a  ****** 
********  fortran  stop  may  be  encountered.  the  amount  of  ****** 
********  work  space  allocated  in  RWKSP  and  IWKIN  must  be  ****** 
********  increased 

COMMON/WORKS P/RWKSP 
REAL  RWKSP  (900000) 

CALL  IWKIN  (900000) 

********  ensure  the  proper  data  files  are  being  read  ********** 
OPEN(7,FILE-'HEAT.OUT' , STATUS- ' NEW' ) 

OPEN(8,FILE-'TESTl,OUT' , STATUS- ' NEW' ) 

0PEN(9,FILE-'TEST2.0UT' , STATUS- ' NEW' ) 

***  reads  given  data  and  creates  mesh 
CALL  INPUT2 
CALL  GRIDPLUS 

***  creates  biga  matrix  and  bigf  vector 
CALL  HEATMAT 

***  replaces  Galerkln  equation  with  direchlet  bdry  temperatures 

***  modifies  equations  so  that  heat  flux  at  boundaries  of  material 

***  interfaces  are  equal 
CALL  HEATMOD 

***  solves  simultaneous  system  of  equations 
CALL  LSARG( NSNP , BIGA, NSNP , BIGF , 1 , TEMP ) 


***  creates  output  file 
999  CALL  OUTPUT 


END 


81 


SUBROUTINE  INPUT2 


*********************************************************************** 
****  Reads  in  data  to  construct  thermal  grid 

INCLUDE  'HEATCOMN.FOR' 

OPEN (11, FILE- 'data8.DAT' , STATUS- 'OLD' ) 

OPEN( 21 , FILE-'TEST.OUT' , STATUS- ' NEW' ) 

READ (11,*)  NCOLl,NCOL2,NROW 
READdl,*)  (XPOS(I)  ,I-1,NDIFX) 

READdl,*)  (  YPOS(  I  )  ,  I-1,NDIFY) 

READdl,*)  (XINCOL(I)  ,  I-l ,  NDI FX-1 ) 

READdl,*)  (YINCRROW(I)  ,  I-l ,  NDI  FY-1 ) 

READdl,*)  (THERCON(I)  ,RHO(I)  ,  SPHT(  I )  ,  I-l ,  NMATL ) 

READ( 11 , * )  DIVSUB,DIVXA1R,DIVYAIR,DIVXS0L 
READ (11,*)  NBCLOW , NBCUP , NBCVERT , NBCSUB , TAMB 
READdl,*)  BLOW,  HUP,  HVERT,HSUB 
READdl,*)  GEN,QTPR 

READ (11,*)  NBDUP , NBDRDEV , NBDUSUB , NBDRI T , NBDLOW 
READ (11,*)  TUP , TRTDEV , TUP SUB , TRIGHT , TLOW 

END 

SUBROUTINE  GRIDPLUS 

******#************************ilr*i^***************Jt************ 

****  overlays  a  grid  on  a  square  or  rectangle  with  up  to  nine 
****  different  mesh  densities  in  X  and  Y  directions  ,  sets  up 
****  geometric  and  material  properties  for  the  Galerkin  eq. 

INCLUDE  'HEATCOMN.FOR' 

C  open( 21 , file- ' test .out ' , status-' new' ) 

DATA  KEL , XVALUE , YVALUE  /NEL* 0,0,0/ 

RBMAXl  -  NCOLl  +  1 
RBMAX2  -  NCOL2  +  1 
RTMAX  -  NROW  +  1 
XCOUNT  -  1 
NN  -  1 
RBMAX  -  RBMAXl 


****  assigns  X  and  Y  coordinates  to  system  Global  nodal  points 

DO  100  I  -  1,NSNP 

X ( I )  -  XVALUE 
Y(I)  -  YVALUE 

I F ( XVALUE . LT . XPOS ( 2 ) . AND . XVALUE . GE . XPOS (1 ) ) THEN 
XINCCOL  -  XINCOL(l) 

ELSEI F( XVALUE.lt. XPOS ( 3 ) . AND . XVALUE . GE . XPOS ( 2 ) )THEN 
XINCCOL  -  XINCOL(2) 

ELSEI F ( XVALUE . LT . XPOS ( 4 ) .AND . XVALUE . GE . XPOS ( 3 ) ) THEN 
XINCCOL  -  XINCOL(3) 

ELSEI F ( XVALUE . LT . XPOS ( 5 ) . AND . XVALUE . GE . XPOS ( 4 ) ) THEN 
XINCCOL  -  XINCOL(4) 

ELSEIF( XVALUE. LT.XPOS( 6) . AND . XVALUE . GE . XPOS ( 5 ) )THEN 
XINCCOL  -  XINCOL(5) 
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ELSEI F ( XVALUE . LT . XPOS ( 9 ) . AND . XVALUE . GE . XPOS ( 8 ) ) THEN 
XINCCOL  -  XINC0L(8) 

ELSE 

XINCCOL  -  XINCOL(9) 

ENDIF 

XVALUE  -  XVALUE  +  XINCCOL 

I F ( YVALUE . LT . YPOS ( 2 ) . AND . YVALUE . GE . YPOS ( 1 ) ) THEN 
YINCROW  -  YINCRROW(l) 

ELS  E I F ( YVALUE . LT . YPOS ( 3 ) . AND . YVALUE . GE . YPOS ( 2 ) ) THEN 
YINCROW  -  YINCRROW(2) 

ELSEIF( YVALUE. LT. YPOS ( 4 ) .AND . YVALUE . GE . YPOS ( 3 ) )THEN 
YINCROW  -  YINCRROW(3) 

ELSEI F ( YVALUE . LT . YPOS ( 5 ) . AND . YVALUE . GE . YPOS ( 4 ) ) THEN 
YINCROW  -  YINCRROW(4) 

ELSEI F ( YVALUE . LT . YPOS ( 6 ) .AND . YVALUE . GE . YPOS ( 5 ) ) THEN 
YINCROW  -  YINCRROW(5) 

ELS  El F ( YVALUE . LT . YPOS ( 7 ) . AND . YVALUE . GE . YPOS ( 6 ) ) THEN 
YINCROW  -  YINCRROW(6) 


ELS  E I F ( YVALUE . LT . YPOS ( 8 ) . AND . YVALUE . GE . YPOS ( 7 ) ) THEN 
YINCROW  -  YINCRROW(7) 

ELSEI F ( YVALUE . LT . YPOS ( 9 ) .AND . YVALUE . GE . YPOS ( 8 ) ) THEN 
YINCROW  -  YINCRROW(8) 

ELSE 

YINCROW  -  YINCRROW(9) 

ENDIF 

****  begins  numbering  process  of  assembly  after  substrate 
IFd.GT.  (  (NCOLl+l)*(DIVSUB+l)  )  )THEN 
RBNAX  -  RBMAX2 
ENDIF 

I F ( XCOUNT . GE . RBMAX ) THEN 
XVALUE  -  0.0 

YVALUE  -  YVALUE  +  YINCROW 
XCOUNT  -  0.0 
ENDIF 

XCOUNT  -  XCOUNT  +  1 
100  CONTINUE 

***************************************************************** 
****  creates  correspondence  table  for  Local-Global  nodal  points 

DO  200  IROW  -  1,DIVSUB  +  1 
KK  -  (IROW-l)*(RBMAXl ) 

DO  210  ICOL  -  l,NCOLl 

lEL  -  ( (IROW-l)*(2*NCOLl) )+(2*ICOL)-l 
JEL  -  lEL  +  1 
ICOR(IEL,l)  -  ICOL  +  KK 
ICOR(IEL,2)  -  ICOL  +  NCOLl  +  KK  +  2 
ICOR(IEL,3)  -  ICOL  +  NCOLl  +  KK  +  1 
ICOR(JEL,l)  -  ICOL  +  KK 
ICOR(JEL,2)  -  ICOL  +  KK  +  1 
ICOR(JEL,3)  -  ICOL  +  NCOLl  +  KK  +  2 
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210  CONTINUE 

200  CONTINUE 


****  correspondence  starting  two  rows  above  substrate 
KK  -  (RBMAX1*(DIVSUB+1) ) 
lEL  -  (NCOLl*DIVSUB*2)  +  (NCOL2*2)  +  1 
DO  220  IROW  -  l,NROW-(DIVSUB+l) 

DO  230  ICOL  -  l,NCOL2 


JEL  -  lEL  +  1 
ICOR(IEL,l)  -  ICOL 
ICOR(IEL,2)  -  ICOL 
ICOR(IEL,3)  -  ICOL 
ICOR(JEL,l)  -  ICOL 
ICOR(JEL,2)  -  ICOL 
ICOR(JEL,3)  -  ICOL 
lEL  -  lEL  -t-  2 


+  KK 

+  NCOL2  +  KK  +  2 
+  NCOL2  +  KK  +  1 
+  KK 

+  KK  +  1 

+  NCOL2  +  KK  +  2 


230  CONTINUE 

KK  -  KK  RBMAX2 
220  CONTINUE 


*************************************************************** 

****  assigns  material  properties  to  elements 

NSUBEL  -  NCOLl  *  DIVSUB  *  2 
NAIREL  -  2  *  DIVXAIR 

NARSOEL  -  2  *  DIVYAIR  *  NCOL2 

LCOUNT  -  1 

NCOUNT  -  1 

MCOUNT  -  1 


DO  300  I  -  1,NEL-1,2 
IF( I .LT. NSUBEL) THEN 
****  substrate  (single  material)  **** 

IMAT(I)  -  1 
IMAT(I+1)  -  1 

****  substrate  (multiple  material)  **** 

ELSEIF( I .GT.NSUBEL.AND. I .LT. ( NARSOEL+NSUBEL ) )THEN 
I F ( NCOUNT . LE . DIVXAIR ) THEN 
****  ai r  **** 

IMAT(I)  -  2 
IMAT(I+1)  -  2 

ELSEI F ( NCOUNT . LE . DI VXAIR+DIVXSOL ) THEN 

****  solder  **** 

IMAT(I)  -  3 
IMAT(I+1)  -  3 
****  coating/cover  **** 

ELSE 

IMAT(I)  -  5 
IMAT(I+1)  -  5 
ENDIF 

NCOUNT  -  NCOUNT  +  1 
I F ( NCOUNT . GT . NCOL2 ) THEN 
NCOUNT  -  1 
ENDIF 

ELSEI F ( I . GT . NARSOEL+NSUBEL .AND . 

:  I .LT.NEL-(NCOL2*2) )THEN 

IF(MCOUNT.LE.NCOL2-l )THEN 
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****  dovico  **** 

IMAT(I)  -  4 
IMAT(I+1)  -  4 
****  coating/cover  **** 

ELSE 

IMAT(I)  -  5 
IMAT(I+1)  -  5 
ENDIF 

MCOUNT  -  MCOUNT  +  1 
I F ( MCOUNT . GT . NCOL2 ) THEN 
MCOUNT  -  1 
ENDIF 
ELSE 

****  coating/cover  **** 

IMAT(I)  -  5 
IMAT(I+1)  »  5 
ENDIF 
"300  CONTINUE 

DO  310  I  -  1,NEL 

II  -  IMAT(I) 

KEL(I)  -  THERCON(II) 
RHOEL(I)  -  RHO(II) 


SPHTEL(I)  -  SPHT(II) 

I F ( I . GT . ( NSUBEL+NARSOEL ) ) THEN 
QDOTEL(I)  -  QTPR 
ENDIF 
310  CONTINUE 

*************************************************************** 

****  assigns  boundary  conditions  to  elements  at  the  convective 

****  (CAUCHY)  boundaries. 

DO  400  I  -  2,2*NC0L1,2 

***  assigns  Cauchy  boundary  to  lower  side  of  rectangle 
I F ( NBCLOW . EQ , 1 ) THEN 
NBDRY(I)  -  1 
ENDIF 
400  CONTINUE 

DO  410  I  -  (NEL-(2*NCOL2) )+l,NEL-l,2 

***  assigns  Cauchy  boundary  to  upper  side  of  rectangle 
IF(NBCUP.EQ.1)THEN 
NBDRY^I)  -  2 
ENDIF 
410  CONTINUE 

DO  420  I  -  NSUBEL-( {NCOLl-NCOL2)*2)+l,NSUBEL-l,2 

***  assigns  Cauchy  boundary  to  upper  tip  of  substrate 
I F ( NBCSUB . EQ . 1 ) THEN 
NBDRYd)  -  4 
ENDIF 
420  CONTINUE 


DO  430  I  -  NSUBEL+{NCOL2*2) ,NEL,NCOL2*2 

***  assigns  Cauchy  boundary  to  vertical  side  of  rectangle 
I F ( NBCVERT . EQ . 1 ) THEN 
NBDRY(I)  -  3 
ENDIF 
430  CONTINUE 

************************************************************** 
****  assigns  Direchlet  boundary  to  system  Global  nodal  points 

ft  it  "k  It 

****  assigns  Direchlet  boundary  to  lower  side  of  substrate 

DO  500  I  -  1,RBMAX1 

I F ( NBDLOW . EQ . 1 ) THEN 
NDIRECH(l)  -  1 
BDRYVAL(I)  -  TLOW 
ENDIF 

500  CONTINUE 

****  assigns  Direchlet  boundary  to  right  hand  side  of  substrate 

DO  510  I  -  RBMAX1*2,RBMAX1*(DIVSUB+1) ,RBMAX1 
I F ( NBDRIT . EQ . 1 ) THEN 
NDIRECH(I)  -  1 


BDRYVAL(I)  -  TRIGHT 
ENDIF 
510  CONTINUE 

****  assigns  Direchlet  boundary  to  upper  side  of  device 

DO  520  I  -  NSNP,NSNP-NCOL2,-l 
I F ( NBDUP . EQ . 1 ) THEN 
NDIRECH(I)  -  1 
BDRYVAL(I)  -  TUP 
ENDIF 
520  CONTINUE 

****  assigns  Direchlet  boundary  to  right  side  of  device 
III  -  NSNP-(RBMAX2*(NROW-DIVSUB-l) ) 

DO  530  I  -  NSNP-RBMAX2,III,-RBMAX2 
I F ( NBDRDEV . EQ . 1 ) THEN 
NDIRECH(I)  -  1 
BDRYVAL(I)  -  TRTDEV 
ENDIF 
530  CONTINUE 

****  assigns  Direchlet  boundary  to  upper  side  of  substrate 

DO  540  I  -  (RBMAXl*DIVSUB)+RBMAX2,RBMAXl*(DIVSUB-»-l  ) 

I F ( NBDUSUB , EQ . 1 ) THEN 
NDIRECH(I)  -  1 
BDRYVAL(I)  -  TUPSUB 
ENDIF 
540  CONTINUE 

END 
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SUBROUTINE  HEATNAT 

************************************************************* 
****  calculates  the  BIGA  matrix  and  BIGF  vector  that  are  the 
****  system  of  simultaneous  equations  to  be  solved;  accounts 
****  for  the  main  operator(Del“2) ,  Cauchy  boundary,  and  the 
****  forcing  function  to  Poisson's  equation 

INCLUDE  ' HEATCOMN . FOR ' 

DATA  BIGA, BIGF  /  NDIM*0,NSNP*0  / 

DATA  bb21,bb22,bb23,bb33  /O . 0 , 0 . 0 , 0 . 0 , 0 . 0/ 

DATA  cauchf 1 , cauchf 2 , f f 1  /O . 0 , 0 . 0 , 0 . 0/ 

DO  10  lEL  -  1,NEL 

11  -  ICOR(IEL,l) 

12  -  ICOR(IEL,2) 

13  -  ICOR(IEL,3) 

B1  -  Y(I2)  -  Y(I3) 

B2  -  Y(I3)  -  Y(I1) 

B3  -  Y(I1)  -  Y(I2) 

Cl  -  X( 13 )  -  X( 12 ) 

C2  -  X(I1)  -  X(13) 

C3  -  X(I2)  -  X(I1) 

AREA(IEL)  -  (ABS(C1*B2-C2*B1) )/2. 
write(8,*)  iel , area( iel ) 

aall  -  KEL( IEL)*(B1**2  +  Cl**2 )/( -4 . *AREA( IEL) ) 
aal2  -  KEL( IEL)*(B1*B2  +  Cl *C2 )/( -4 . *AREA( IEL ) ) 

aal3  =•  KEL(IEL)*(B1*B3  +  C1*C3  )/( -4  .  *AREA(IEL )  ) 

aa21  -  KEL(IEL)*(B2*B1  +  C2*C1 )/( -4 . *AREA( IEL ) ) 

aa22  -  KEL( IEL ) * ( B2**2  +  C2**2 )/( -4 . *AREA( IEL ) ) 
aa23  -  KEL( IEL) *(B2*B3  +  C2*C3 )/( -4 . *AREA( IEL ) ) 

aa31  -  KEL( IEL) *(B3*B1  +  C3*Cl )/( -4 . *AREA( IEL ) ) 

aa32  -  KEL(IEL)*(B3*B2  +  C3*C2 )/( -4 . *AREA( IEL ) ) 

aa33  -  KEL( IEL ) * ( B3**2  +  C3**2 )/( -4 . *AREA( IEL ) ) 


************************************************************** 
****  calculation  of  Cauchy  boundary  to  BIGA  matrix  and  BIGF 
****  vector 

IF(NBDRy( IEL) .NE. 0 )THEN 

****  boundary  with  local  nodal  points  1  -  2  on  bottom  side 

IF{NBDRY( IEL) . EQ. 1 )THEN 

bb22  -  (HL0W*ABS(C3) )/-3.0 

bb21  -  (HLOW*ABS(C3) )/-6.0 

cauchf 1  -  (TAMB*HLOW*ABS(C3) )/-2.0 
****  boundary  with  local  nodal  points  2  -  3  on  upper  side 

ELSEIF(NBDRY( IEL) .EQ.2)THEN 

bb33  -  (HUP*ABS(C1) )/-3.0 

bb23  -  {HUP*ABS(C1) )/-6.0 

cauchf 2  -  (TAMB*HUP*ABS(C1 ) )/-2.0 

****  boundary  with  local  nodal  points  2  -  3  on  right  side 

ELSEIF(NBDRY( IEL) .EQ.3)THEN 

bb33  -  (HVERT*ABS(B1) )/-3.0 

bb23  -  (HVERT*ABS(B1) )/-6.0 
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cauchf2  -  (TAMB*HVERT*ABS(B1) )/-2.0 
ELSE 

bb33  -  (HSUB*ABS(C1) )/-3.0 
bb23  -  (HSUB*ABS(C1) )/-6.0 
cauchf2  -  (TAMB*HSUB*ABS(C1) )/-2.0 
END  IF 
END  IP 

*********************************************************** 

****  the  contribution  of  the  Cauchy  boundary  is  symmetric 
****  bbll-bb22  bbl2-bb21  bb22-bb33  bb23-bb32 

BIGA(I1,I1)  -  BIGA(I1,I1)  +  aall  +  bb22 

BIGA(I1,I2)  -  BIGA(I1,I2)  +  aal2  +  bb21 

BIGA(I1,I3)  -  BIGA(I1,I3)  +  aal3 

BIGA(I2,I1)  -  BIGA(I2,I1)  +  aa21  +  bb21 

BIGA{I2,I2)  -  BIGA(I2,I2)  +  aa22  +  bb22  +  bb33 

BIGA(I2,I3)  -  BIGA(I2,I3)  +  aa23  +  bb23 

BIGA(I3,I1)  -  BIGA(I3,I1)  +  aa31 

BIGA(I3,I2)  -  BIGA(I3,I2)  +  aa32  +  bb23 

BIGA(I3,I3)  -  BIGA(I3,I3)  +  aa33  +  bb33 

bb22  -  0.0 
bb21  -  0.0 
bb23  -  0.0 
bb33  -  0.0 

****  used  for  Poisson's  equation  (internal  heat  generatic:* ; 

IF(GEN.EQ.1.0)THEN 

IF( IMAT( lEL) .EQ. 4 )THEN 

ffl  =  -AREA(IEL)*QTPREL{IEL)/3.0 
ENDIF 

ENDIF 

****  the  two  terms  that  comprise  the  BIGF  matrix  are  from  the 
****  Cauchy  boundary  condition  and  the  source  term  of  the  chip 
****  a  lumped  approximation  is  used  for  both  cases;  cauchfl  is 
****  for  a  1-2  boundary  edge  cauchf2  is  for  a  2-3  boundary  edge 

BIGF(ll)  -  BIGF(ll)  +  ffl  +  cauchfl 

BIGF(I2)  -  BIGF(I2)  +  ffl  +  cauchfl  +  cauchf2 

BIGF(I3)  -  BIGF(I3)  +  ffl  +  cauchf2 

cauchfl  -  0.0 

cauchf2  -  0.0 

ffl  -  0.0 

10  CONTINUE 
END 

SUBROUTINE  HEATMOD 

****************************************************************** 
****  modifies  BIGA  matrix  and  BIGF  vector  if  a  system  Global  nodal 
****  point  is  a  Direchlet  boundary  ie.  replaces  the  Galerkin  eq. 
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INCLUDE  ' HBATCOMN . FOR 

****  modifies  system  of  equations  for  Direchlet  boundary  points 
DO  10  I  -  1,NSNP 

IF  (NDIRECH(I)  .EQ.  1)THEN 
BIGA(I,I)  -  1.0 

DO  20  J  -  1,1-1 
BIGA(I,J)  -  0.0 
20  CONTINUE 

DO  30  JJ  »  I+1,NSNP 
BIGA(I,JJ)  -  0.0 
30  CONTINUE 

BIGF(I)  «  BDRYVAL(I) 

ENDIF 

10  CONTINUE 
END 

SUBROUTINE.  OUTPUT 

******************************************************************** 
****  creates  output  file  of  temperature  distribution  and  a  data  file 
****  for  graphing  the  temperature  contours 

INCLUDE  'HEATCOMN.FOR- 

OPEN( 10, FILE-'HEAT.DAT' , STATUS-'NEW' ) 
write(7,*'  '  Temperatures  through  solder' 

wrlte(7,*)  '8', temp{ 8 ) , ' 23 ' , temp( 23 ) , ' 38 ' , temp( 38 ) , ' 53 ' , temp( 53 ) 
write ( 7 , * )  ' 68 ' , temp ( 68 ) , ' 83 ' , temp( 83 ) , ' 93 ' , temp( 93 ) , 

!  ' 103 ' , temp( 103 ) 

write (7,*)  '113' ,temp(113) , ' 123 ' , temp( 123 ) , ' 133 ' ,temp(133) 

****  energy  balance 

qconv  -  0.0 

do  200  i  -  126,134 

qconv-qconv+( hup* ( ( ( tamp( I )+temp( I+l ) )/2 )-tamb) * (X(I+l)-x(I))) 

200  continue 

do  210  i  -  1,14 

qconv-qconv+( blow* ( ( ( temp( I )+temp( I+l ) )/2 )-tamb) * (X( I+l )-x( I ) ) ) 
210  continue 

do  220  i  -  70,74 

qconv«-qconv+  ( hsub*  (  ( ( terap(  I  )+temp(  I  +  l )  )/2  )-tamb )  *  (X(  I+l  )-x(  I )  )  ) 

.  220  continue 

do  230  i  -  85,125,10 

qconv-qconv+(hvert*( ( ( temp( I )+temp( I+IO) )/2)-tamb)* 
s  (y(I+10)-y( I) ) ) 

230  continue 

qconv-qconv+( hvert* ( ( ( temp( 70 ) 4  terop( 85 ) )/2 )-tamb) * (y( 85)-y ( 70 ) ) ) 
WRITE( 7 , * ) 

write(7,*)  '  Heat  loss  to  convection  -  ', qconv, 'W/M' 

POWER  -  QPRT*0. 0025*0. 002 

WRITE(7,*)  'The  power  output  of  the  device  -', power , 'W/M' 

END 
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APPENDIX  B. 


PROGRAM  FILTER 

REAL  X(135),Y(135),T(135) ,TEMPT(200) ,TEMPB(200) 

OPEN{12,FILE-'bh3.out' , STATUS-'OLD' ) 

OPEN( 13, FILE- 'bh3f.DAT' , STATUS- ' NEW' ) 

READ{12,*)  (X( I) ,Y(I) ,T(I) ,1-1,135) 

****  calculates  intermediate  temperatures  ******* 

TA  -  (T(61)-T(46) )*.75  +  T(46) 

TB  -  (T(61)+T(46) )/2.0 

TC  -  (T(62)-T( 47) )*.75  +  T(47) 

TD  -  (T(66)+T( 47) )/2.0 

TE  -  (T(63)-T(48) )*.75  +  T(48) 

TF  -  (T(63)+T( 48) )/2.0 

TG  -  (T(64)-T(49) )*.75  +  T(49) 

TH  -  (T(64)+T(49) )/2.0 

TI  -  (T(65)-T(64) )*.75  +  T(64) 

TL  -  (T( 50)-T( 49) )*.75  +  T(49) 

TJ  -  (TI-TL)*.75  +  TL 

TK  -  (TI+TL)/2,0 

TM  -  (T(65)-T(50) )*.75  +  T(50) 

TN  -  (T(65)+T(50) )/2.0 
TO  -  (T(65)+T(66) )/2.0 
TR  -  (T(51)-T(50) )*.75  +  T(50) 

TP  -  (TO-TR)*,75  +  TR 

TQ  -  (TO+TR)/2,0 

TS  -  (T(65)+T(66) )/2,0 

TV  -  (T(50)+T(51) )/2.0 

TT  -  (TS-TV)*.75  +  TV 

TU  -  (TS+TV)/2.0 

TW  -  {T(66)-T(51) )*.75  +  T(51) 

TX  -  (T(66)+T(51) )/2.0 

TY  -  (T(67)-T(52) )*.75  +  T(52) 

TZ  -  (T(67)+T(52) )/2.0 

TAA  -  (T(68)-T(53) )*.75  +  T(53) 

TAB  -  (T{68)+T( 53) )/2.0 

TAG  -  (T(69)-T(54) )*.75  +  T(54) 

TAD  -  (T(69)+T(54) )/2.0 

TAE  -  (T(69)+T(70) )/2.0 

TAH  -  (T(54)+T(55) )/2.0 

TAF  -  (TAE-TAH)*.75  +  TAH 

TAG  -  (TAE+TAH)/2.0 

TAI  -  (T(69)-T(70) )*.75  +  T(70) 

TAL  -  (T(54)-T(55) )*.75  +  T(54) 

TAJ  -  (TAI-TAL)*.75  +  TAL 

TAX  -  (TAI+TAL)/2.0 

TAM  -  (T(70)-T(55) )*.75  +  T(55) 

TAN  -  (T(70)+T(55) )/2.0 

TAO  -  (T(70)-T(71) )*.75  +  T(71) 
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TAR  -  (T( 55)-T{ 56) )*.75  +  T(55) 

TAP  -  (TAO-TAR)*.75  +  TAR 

TAQ  -  (TAO+TAR)/2.0 

TAS  -  {T(71)-T(56) )*.75  +  T(56) 

TAT  -  (T(71)+T(56) )/2.0 

TAU  -  (T(75)-T(60) )*.75  +  T(60) 

TAV  -  (T(75)+T(60) )/2.0 

TAW  -  (T(35)-T(34) )*.75  +  T(34) 

TAX  -  (T(36)-T{35) )*.75  +  T(35) 

TAY  -  (T(36)+T(35))/2.0 

TAZ  -  (T( 39)+T(40) )/2.0 

TEA  -  (T(39)-T(40) )*.75  +  T(40) 
TBB  -  (T(40)-T(41) )*.75  +  T(41) 
TBC  -  (T(5)-T(4) )*.75  +  T(4) 

TBD  -  (T(6)-T(5) )*.75  +  T(5) 

TBE  -  (T(6)+T(5) )/2.0 

TBP  -  (T(9)+T(10) )/2.0 

TBG  -  (T(9)-T(10) )*.75  +  T(10) 

TBH  -  (T(lO)-T(ll) )*.75  +  T(ll) 


****  calculates  average  temperatures  for  rectanqular 
****  weld. for  (chip  and  solder) 


element  FEM 


program 


TEMPSl  -  (T(95)+T(85)+T(96)+T(86) )/4.0 
DO  100  I  -  1,115 

TEMPT(I)  -  TEMPSl 
TEMPB(I)  -  TEMPSl 
100  CONTINUE 


!!!!  <^3lculates  average  temperatures  for  rectangular 
****  weld. for  (substrate) 


element 


FEM 


program 


TEMPT(116)  «  (T(61)+T(61) )/2.0 

TEMPB(116)  -  (TA+TC)/2.0 


TEMPT(117)  -  (T(62)+T(63) )/2.0 

TEMPB(117)  «  (TC+TD)/2.0 


TEMPT(118)  «  (T(63)+T(64) )/2.0 

TEMPB(118)  «  (TE+TG)/2.0 


TEMPT(119)  »  (T(64)+TI)/2.0 

TEMPB(119)  -  (TG+TJ)/2.0 

TEMPT(120)  -  (T(65)+TI )/2.0 

TEMPB(120)  -  (TJ+TM)/2.0 

TEMPT(121)  -  (T(65)+TO)/2.0 

TEMPB(121)  -  (TM+TP)/2.0 

TEMPT(122)  -  (TO+TS)/2.0 

TEMPB(122)  -  (TP+TT)/2.0 

TEMPT(123)  -  (T(66)+TS)/2.0 

TEMPB(123)  -  (TT+TW)/2.0 


TEMPT(124)  -  (T(66)+T(67) )/2.0 

TEMPB(124)  -  (TW+TY)/2.0 
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r 


TEMPT(125) 

-  (T(67)+T(68) )/2.0 

TEMPB( 125) 

-  (TY+TAA)/2.0 

TEMPT(126) 

-  {T(68)+T(69) )/2.0 

TEMPB(126) 

-  (TAA+TAC)/2.0 

TEMPT(127) 

-  (T(69)+TAE)/2.0 

TEMPB(127) 

-  (TAC+TAF)/2.0 

TEMPT(128) 

-  (TAE+TAI )/2.0 

TEMPB(128) 

-  (TAF+TAJ)/2.0 

TEMPT( 129) 

-  (T(70)+TAI)/2.0 

TEMPBU29) 

-  (TAJ+TAM)/2.0 

TEMPT(130) 

-  (T(70)+TAO)/2.0 

TEMPB(130) 

-  (TAM+TAP)/2.0 

TEMPT( 131) 

-  (T( 71 )+TAO)/2 .0 

TEMPB( 131 ) 

-  (TAP+TAS)/2.0 

TEMPT( 132) 

-  (T(71)+T(75) )/2.0 

TEMPB(132) 

-  (TAS+TAU)/2.0 

TEMPT( 133 ) 

-  (TA+TC)/2.0 

TEMPB( 133  ) 

-  (TB+TD)/2.0 

TEMPT( 134 ) 

-  (TC+TE)/2.0 

TEMPB(134) 

«  (TD+TF)/2.0 

TEMPT( 135) 

-  (TE+TG)/2.0 

TEMPB(135) 

-  (TF+TH)/2.0 

TEMPT(136) 

-  (TG+TJ)/2.0 

TEMPB( 136) 

-  (TH+TK)/2.0 

TEMPT(137) 

-  (TJ+TM)/2.0 

TEMPB( 137 ) 

-  (TK+TN)/2.0 

TEMPT(138) 

-  (TM+TP)/2.0 

TEMPB(138) 

-  (TN+TQ)/2.0 

TEMPT(139) 

-  (TP+TT)/2.0 

TEMPB(139) 

-  (TQ+TU)/2.0 

TEMPT( 140) 

-  (TT+TW)/2.0 

TEMPB( 140) 

-  (TU+TX)/2.0 

TEMPT(141) 

-  (TW+TY)/2.0 

TEMPB( 141 ) 

-  (TX+TZ)/2.0 

TEMPT( 142 ) 

-  (TY+TAA)/2.0 

TEMPB(142) 

-  (TZ+TAB)/2.0 

TEMPT(143) 

-  (TAA+TAC)/2.0 

TEMPB(143) 

-  (TAB+TAD)/2.0 

TEMPT(144) 

TEnPBa44) 


{TAC+TAF)/2.0 

(TAD+TAG)/2.0 
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TEMPT( 145) 
TEMPB( 145) 


(TAF+TAJ)/2.0 

(TAG+TAK)/2.0 


TEMPT(146)  -  (TAJ+TAM)/2.0 

TEMPB(146)  -  (TAK+TAN)/2.0 

TEMPT(147)  -  (TAM+TAP)/2.0 

TEMPB(147)  -  (TAN+TAQ)/2.0 

TEMPT(148)  -  (TAP+TAS)/2.0 

TEMPB(148)  -  (TAQ+TAT)/2.0 

TEMPT(149)  -  (TAS+TAU)/2 .0 

TEMPB(149)  -  (TAT+TAV)/2.0 

TEMPT(150)  -  (TB+TD)/2.0 

TEMPB(150)  -  (T(46)+T(47) )/2.0 

TEMPTdSl)  =  (TD+TF)/2.0 

TEMPB(151)  =■  (T(47)+T(48)  )/2.0 

TEMPT(152)  -  (TF+TH)/2.0 

TEMPB(152)  -  (T(48)+T(49) )/2.0 

TEMPT(153)  -  (TH+TK)/2.0 

TEHPB(153)  -  (TL+T( 49) )/2.0 

TEMPT(154)  -  (TK+TL)/2.0 

TEMPB(154)  -  (TN+T( 50) )/2.0 

TEMPT(155)  -  (TN+TQ)/2.0 

TEMPB(155)  -  <TR+T(50) )/2.0 

TEMPT(156)  «  (TQ+TU)/2.0 

TEMPB(156)  «  (TR+TV)/2.0 

TEMPT(157)  -  (TU+TX)/2.0 

TEMPB(157)  -  (TV+T(51) )/2.0 

TEMPT(158)  -  (TX+TZ)/2.0 

TEMPB(158)  -  (T(51)+T(52) )/2.0 

TEMPT (159)  -  (TZ+TAB)/2.0 

TEMPB(159)  -  (T(52)+T(53) )/2.0 

TEMPT(160)  -  (TAB+TAD)/2.0 

TEMPB(160)  -  (T(53)+T(54) )/2.0 

TEMPT(161)  -  (TAD+TAG)/2.0 

TEMPBaei)  -  (T(54)+TAH)/2,0 

TEMPT(162)  -  (TAG+TAK)/2.0 

TEMPB(162)  -  (TAH+TAL)/2.0 

TEMPT(163)  -  (TAK+TAN)/2.0 

TEMPB(163)  -  (TAL+T( 55) )/2.0 

TEMPT(164)  -  (TAN+TAQ)/2.0 

TEMPB(164)  -  {TAR+T{ 55) )/2.0 
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TEMPT(165)  -  (TAQ+TAT)/2.0 

TEMPB(165)  -  (TAR+T(56) )/2.0 

TEMPT(166)  -  (TAT+TAV)/2.0 

TEMPB(166)  -  (T( 56)+T(60) )/2.0 

TEMPT{167)  -  (T(46)+T(47) )/2.0 

TEMPB(167)  -  (T( 31)+T( 32) )/2.0 

TEMPT(160)  -  (T(47)+T( 48) )/2.0 

TEMPB(168)  -  {T( 32)+T( 33) )/2.0 

TEMPT(169)  -  (T( 48)+T( 49) )/2.0 

TEMPB(169)  -  (T(33)+T(34) )/2.0 

TEMPT(170)  -  (T(49)+TL)/2.0 

TEMPB(170)  »  (T( 34)+TAW)/2.0 

TEMPT(171)  -  (T(50)+TL)/2.0 

TEMPB(171)  -  (T( 35)+TAW)/2.0 

TEMPT(172)  -  (T(50)+TR)/2.0 

TEMPB(172)  -  (T( 35)+TAX)/2.0 

TEMPT(173)  »  (TR+TV)/2.0 

TEMPB(173)  -  (TAX+TAY)/2.0 

TEMPT(174)  -  (T(51)+TV)/2.0 

TEMPB(174)  -  (T(36)+TAY)/2,0 

TEMPT(175)  -  (T(51)+T(52) )/2.0 

TEMPB(175)  -  (T(36)+T(37) )/2.0 

TEMPT(176)  -  {T( 52)+T( 53) )/2.0 

TEMPB(176)  -  (T(37)+T(38))/2.0 

TEMPT(177)  -  (T(53)+T(54) )/2.0 

TEMPB(177)  -  (T( 38)+T( 39) )/2.0 

TEMPT(178)  -  (T( 54 )+TAH)/2.0 

TEMPB(178)  -  {T(39)+TAZ)/2.0 

TEMPT(179)  -  (TAH+TAL)/2.0 

TEMPB(179)  -  (TAZ+TBA)/2.0 

TEMPT(180)  -  (T(55)+TAL)/2.0 

TEMPB(180)  -  (T(40)+TBA)/2.0 

TEMPT(181)  -  (T(55)+TAR)/2.0 

TEMPB(181)  -  (T(40)+TBB)/2.0 

TEMPT(182)  -  (T(56)+TAR)/2.0 

TEMPB(182)  -  (T( 41 )+TBB)/2 .0 

TEMPT(183)  -  (T(56)+T(60) )/2.0 

TEMPB(183)  -  (T(41 )+T( 45) )/2.0 

TEMPT(184)  -  (T(31)+T(32) )/2.0 

TEMPB(184)  -  (T(l)+T(2) )/2.0 
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TEMPT(185) 
TEMPB( 185) 

-  (T(32)+T(33) )/2.0 

-  (T(2)+T(3) )/2.0 

TEMPT(186) 
TEMPB( 186 ) 

-  (T(33)+T(34) )/2.0 

-  (T(3)+T(4) )/2.0 

TEMPT(187) 
TEMPB{ 187 ) 

-  {T(34)+TAW)/2.0 

-  (T(4)+TBC)/2.0 

TEMPT(188) 

TEMPB(188) 

-  (T(35)+TAW)/2.0 

-  (T( 5)+TBC)/2.0 

TEMPT( 189 ) 
TEMPB(189) 

-  (T( 35)+TAX)/2.0 

-  (T( 5)+TBD)/2.0 

TEMPT( 190) 
TEMPB(190) 

-  (TAX+TAY)/2.0 

-  (TBD+TBE)/2.0 

TEMPT( 191 ) 
TEMPB( 191 ) 

-  (T( 36)+TAY)/2.0 

-  (T(6)+TBE)/2.0 

TEMPT( 192) 
TEMPB( 192) 

-  (T(36)+T(37) )/2.0 

-  (T(6)+T(7) )/2.0 

TEMPT(193) 

TEMPB(193) 

-  (T(37)+T(38) )/2.0 

-  (T(7)+T(8) )/2.0 

TEMPT( 194 ) 
TEMPB( 194) 

-  (T( 38)+T( 39) )/2.0 

-  (T(8)+T(9) )/2.0 

TEMPT( 195) 
TEMPB( 195) 

-  (T( 39)+TAZ)/2.0 

-  (T(9)+TBF)/2.0 

TEMPT(196) 

TEMPB(196) 

-  (TAZ+TBA)/2.0 

-  (TBF+TBG)/2.0 

TEMPT( 197 ) 
TEMPB(197) 

-  (T(40)+TBA)/2.0 

-  (T(10)+TBG)/2.0 

TEMPT( 198) 
TEMPB( 198) 

-  (T(40)+TBB)/2.0 

-  (T(10)+TBH)/2.0 

TEMPT(199) 
TEMPB( 199) 

-  (T(41)+TBB)/2.0 

-  (T(ll)+TBH)/2.0 

TEMPT(200) 

TEMPB(200) 

-  (T(41)+T(45) )/2.0 

-  (T(ll)+T(15) )/2.0 

***  use  if  printout  is  desired 

DO  900  I  -  1,200 

WRITE (13,*)  TEMPT (I),'  ',TEMPB(I) 
900  CONTINUE 
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